%Program rkstep.m for Runge-Kutta's 4th order method
%Model problem is Van der Pol's equation, stored in fvdp.m
 global epsilon      %parameter used in both script and fvdp 
 epsilon=1;          %parameter in VdPs equation
 t0=0;tend=10;       %time interval for the solution
 u0=[1 0]';          %initial value
 N=100;              %number of steps
 h=(tend-t0)/N;      %stepsize
 result=[u0'];       %collecting solution values
 time=[t0];          %collecting time points
 u=u0;t=t0;          %initialization of RK's method
 for k=1:N
     k1=fvdp(t,u);
     k2=fvdp(t+h/2,u+h*k1/2);
     k3=fvdp(t+h/2,u+h*k2/2);
     k4=fvdp(t+h,u+h*k3);
     u=u+h*(k1+2*k2+2*k3+k4)/6;
     t=t+h;
     result=[result;u'];
     time=[time;t];
 end
 plot(time,result(:,1),'o',time,result(:,2),'+')
 title('Van der Pols equation')
 xlabel('time t')
 ylabel('u (o) and du/dt (+)')
 grid


%Right hand side of Van der Pol's equation
 function rhs=fvdp(t,u)
 global epsilon     %parameter in VdPs equation
 rhs=[u(2);
      -epsilon*(u(1)*u(1)-1)*u(2)-u(1)];



%Program showing the use of Matlabs ode45 function
%Model problem is Van der Pol's equation, stored in fvdp.m
 global epsilon  
 epsilon=1;         %parameter in VdPs equation
 t0=0;tend=10;      %time interval for the solution
 u0=[1 0]';         %initial value
 [time,result]=ode45(@fvdp,[t0,tend],u0);
 plot(time,result(:,1),'o',time,result(:,2),'+')
 title('Van der Pols equation')
 xlabel('time t')
 ylabel('u (o) and du/dt (+)')
 grid


%Example 1.2 in chapter 1. The vibration equation
%demo of tolerance parameter settings
 global m c k f0 w
 m=1;c=10;k=1e3;f0=1e-4;w=40;
 t0=0;tend=5;
 u0=[0 0]';
 [time,result]=ode45(@fvib,[t0,tend],u0);    %not enough accuracy
 options=odeset('RelTol',1e-6,'AbsTol',1e-9);%set tolerances
 [t,res]=ode45(@fvib,[t0,tend],u0,options);  %accurate solution
 plot(time,result(:,1),'.',t,res(:,1))
 title('Accurate and bad solution of the vibration equuation')
 xlabel('time t [sec]')
 ylabel('vibration y(t) [m]')