%MATLAB examples - April 23, 2009 (2nd. order eqns via 1st order systems)
%EXAMPLE 1: linear second order ODE
%x''+3x'+2x=sin(t),  x=x(t)  x(0)=1,  x'(0)=0,   t in [0,10]
%equivalent system:  x'=y, y'=-3y-2x+sin(t)
function xprime=system(t,x)
xprime=zeros(2,1);
xprime(1)=x(2);
xprime(2)=-3*x(2)-2*x(1)+sin(t);
%one can also try discontinuous forcing (t>=1)-(t>=2) or
%"ramp up" t.*(t<1)+(t>=1)
%[t,x]=ode45(@system,[0,10],[1;0]));
%plot(t,x)  plots both components
%plot(t,x(:,1)) plots first component
%plot(x(:,1),x(:,2))  phase plane plot
%---------------------------
%EXAMPLE 2: motion under central force proportional to distance
%x''=-x/4, y''=-y/4, x(0)=2, x'(0)=0, y(0)=0, y'(0)=2
%xprime=zeros(4,1);
%xprime(1)=x(2);
%xprime(2)=x(1)/4;
%xprime(3)=x(4);
%xprime(4)=x(3)/4;
%[t,x]=ode45(@system,[0,2],[2;0;0;2])
%-----------------------------
%EXAMPLE 3:motion under inverse-square force
%x''=-100x/(x^2+y^2)^(3/2), y''=-100y/(x^2+y^2)^(3/2),x(0)=10,x'(0)=0,y(0)=0,y'(0)=1
%period: T=2pi*a^(3/2)k^(-1/2)[Kepler], so T=(pi/5)a^(3/2)
%semimajor axis: a=(l^2/k)/(1-e^2)=v^2/(1-e^2), e=|v^2/10-1|, v=y'(0)
%xprime=zeros(4,1);
%xprime(1)=x(2);
%xprime(2)=-100*x(1)/(x(1)^2+x(3)^2)^(3/2);
%xprime(3)=x(4);
%xprime(4)=-100*x(3)/(x(1)^2+x(3)^2)^(3/2);
%[t,x]=ode45(@system,[0,9],[10;0;0;2])
%---------------------------------
%EXAMPLE 4: a nonlinear equation x''=2xx', x=x(t),x(0)=1,x'(0)=2
%xprime=zeros(2,1);
%xprime(1)=x(2);
%xprime(2)=2*x(1)*x(2);
%[t,x]=ode45(@system,[0,0.78],[1;2])
%-----------------------------------------------
%EXAMPLE 5: a conservative autonomous equation x''=8x-4x^3,x(0)=1.5,x'(0)=0
%xprime=zeros(2,1);
%xprime(1)=x(2);
%xprime(2)=8*x(1)-4*x(1)^3;
%[t,x]=ode45(@system,[0,2],[1.5;0])