Your problem is stiff. Let's substitute in symbolic values for y and see what your ODE function will return.
>> vpa(beam(1, [sym('y1'); sym('y2')]))
ans =
y2
- 1583163086609297.0*y1 - 79577963.950060874223709106445312*y2
A small change in either y1 or y2 results in a huge change in dy2/dt and thus a huge change in y2 at the next time step. So the ODE solver needs to take very, very small steps. You can see this if you plot your ODE:
[t,y]=ode45(beam,tspan,ic, odeset('OutputFcn', @odeplot));
Let that run for a minute or two and see that progress is being made, but it is being made very slowly. You may want to zoom in around y = 0, x between 0 and 0.01. Now let's try this same experiment but with a stiffer solver like ode15s:
[t,y]=ode15s(beam,tspan,ic, odeset('OutputFcn', @odeplot));
The plot (and the ODE solver) finished too quickly for me to see the individual time steps being plotted except as a bit of flicker.
