%Solution of Second Order DE
%-------------------------
%d2x/dt2+0.1(dx/dt)^2+0.6x=0
%-------------------------

f1=@(t,x1,x2) x2;                   %dx1/dt=x2                  equation I
f2=@(t,x1,x2) -0.1*x2.^2-0.6*x1;    %dx2/dt=-0.1*x2.^2-0.6*x1   equation II
m=200;                              %number of x and t values
t=linspace(0.0,20,m);
x1=zeros(1,m);
x2=zeros(1,m);
x1(1)=1;                            %initial conditions
x2(1)=0;
h=t(2)-t(1);                        %interval

for n=1:m-1                         %4th order R-K Method
    k11=h*f1(t(n),x1(n),x2(n));     %k values for equation I
    k12=h*f1(t(n)+h/2,x1(n)+k11/2,x2(n));
    k13=h*f1(t(n)+h/2,x1(n)+k12/2,x2(n));
    k14=h*f1(t(n)+h,x1(n)+k13,x2(n));

    x1(n+1)=x1(n)+(1/6)*(k11+2*k12+2*k13+k14);

    k21=h*f2(t(n),x1(n),x2(n));     %k values for equation II
    k22=h*f2(t(n)+h/2,x1(n),x2(n)+k21/2);
    k23=h*f2(t(n)+h/2,x1(n),x2(n)+k22/2);
    k24=h*f2(t(n)+h,x1(n),x2(n)+k23);

    x2(n+1)=x2(n)+(1/6)*(k21+2*k22+2*k23+k24);
end;
subplot(1,2,1)
plot(t,x1,t,x2,'LineWidth',2);                % graphs t vs x t vs dx/dt
legend('t vs x','t vs dx/dt','Location','SouthWest');
subplot(1,2,2);
plot(x1,x2,'LineWidth',2);                    %graph x vs dx/dt
legend('X1 vs X2','Location','SouthWest');
clear all;