RK4 s

%Runge-Kutta to solve differential equations.
%minster
% solve x’ = f(t,x); x(a) = xa
%Demo system is f1 = 2; f2 = y2; f3 = -y3; f4 = (5 – y4)*y4
% y1(0)=2; y2(0)=1; y3(0)=5; y4(0)=1
clear all
close all

%%Initialize
a = 0; %initial time
b = 50; %final time
y1a=4; y2a=0; y3a=0; y4a=0.5; %initial solutions

f1 = inline(‘2′,’t’,’y1′,’y2′,’y3′,’y4′);
f2 = inline(‘y2′,’t’,’y1′,’y2′,’y3′,’y4′);
f3 = inline(‘-y3′,’t’,’y1′,’y2′,’y3′,’y4′);
f4 = inline(‘0.1*y4 + y3*y4′,’t’,’y1′,’y2′,’y3′,’y4′);

%n=100; %Number of time steps
h = 0.25%(b-a)/n; %step size
n=(b-a)/h;
t = (a+h:h:b); %vector of times
y1 = t; y2 = t; y3 = t; y4 = t;
%xrg4 = t; %solution vector. Copy time into it to create it with appropriate length.
%yrg4 = xrg4

%% Begin Runge-Kutta 4th order.

q1 = h*f1(a,y1a,y2a,y3a,y4a);
r1 = h*f2(a,y1a,y2a,y3a,y4a);
v1 = h*f3(a,y1a,y2a,y3a,y4a);
w1 = h*f4(a,y1a,y2a,y3a,y4a);

q2 = h*f1(a + h/2, y1a + q1/2, y2a + r1/2, y3a + v1/2, y4a + w1/2);
r2 = h*f2(a + h/2, y1a + q1/2, y2a + r1/2, y3a + v1/2, y4a + w1/2);
v2 = h*f3(a + h/2, y1a + q1/2, y2a + r1/2, y3a + v1/2, y4a + w1/2);
w2 = h*f4(a + h/2, y1a + q1/2, y2a + r1/2, y3a + v1/2, y4a + w1/2);

q3 = h*f1(a + h/2, y1a + q2/2, y2a + r2/2, y3a + v2/2, y4a + w2/2);
r3 = h*f2(a + h/2, y1a + q2/2, y2a + r2/2, y3a + v2/2, y4a + w2/2);
v3 = h*f3(a + h/2, y1a + q2/2, y2a + r2/2, y3a + v2/2, y4a + w2/2);
w3 = h*f4(a + h/2, y1a + q2/2, y2a + r2/2, y3a + v2/2, y4a + w2/2);

q4 = h*f1(a + h, y1a + q3, y2a + r3, y3a + v3, y4a + w3);
r4 = h*f2(a + h, y1a + q3, y2a + r3, y3a + v3, y4a + w3);
v4 = h*f3(a + h, y1a + q3, y2a + r3, y3a + v3, y4a + w3);
w4 = h*f4(a + h, y1a + q3, y2a + r3, y3a + v3, y4a + w3);

y1(1) = y1a +(1/6)*q1 + (1/3)*q2 + (1/3)*q3 + (1/6)*q4;
y2(1) = y2a + (1/6)*r1 + (1/3)*r2 + (1/3)*r3 + (1/6)*r4;
y3(1) = y3a + (1/6)*v1 + (1/3)*v2 + (1/3)*v3 + (1/6)*v4;
y4(1) = y4a + (1/6)*w1 + (1/3)*w2 + (1/3)*w3 + (1/6)*w4;

for i = 2:n
q1 = h*f1(t(i-1),y1(i-1),y2(i-1),y3(i-1),y4(i-1));
r1 = h*f2(t(i-1),y1(i-1),y2(i-1),y3(i-1),y4(i-1));
v1 = h*f3(t(i-1),y1(i-1),y2(i-1),y3(i-1),y4(i-1));
w1 = h*f4(t(i-1),y1(i-1),y2(i-1),y3(i-1),y4(i-1));

q2 = h*f1(t(i-1) + h/2, y1(i-1) + q1/2, y2(i-1) + r1/2, y3(i-1) + v1/2, y4(i-1) + w1/2);
r2 = h*f2(t(i-1) + h/2, y1(i-1) + q1/2, y2(i-1) + r1/2, y3(i-1) + v1/2, y4(i-1) + w1/2);
v2 = h*f3(t(i-1) + h/2, y1(i-1) + q1/2, y2(i-1) + r1/2, y3(i-1) + v1/2, y4(i-1) + w1/2);
w2 = h*f4(t(i-1) + h/2, y1(i-1) + q1/2, y2(i-1) + r1/2, y3(i-1) + v1/2, y4(i-1) + w1/2);

q3 = h*f1(t(i-1) + h/2, y1(i-1) + q2/2, y2(i-1) + r2/2, y3(i-1) + v2/2, y4(i-1) + w2/2);
r3 = h*f2(t(i-1) + h/2, y1(i-1) + q2/2, y2(i-1) + r2/2, y3(i-1) + v2/2, y4(i-1) + w2/2);
v3 = h*f3(t(i-1) + h/2, y1(i-1) + q2/2, y2(i-1) + r2/2, y3(i-1) + v2/2, y4(i-1) + w2/2);
w3 = h*f4(t(i-1) + h/2, y1(i-1) + q2/2, y2(i-1) + r2/2, y3(i-1) + v2/2, y4(i-1) + w2/2);

q4 = h*f1(t(i-1) + h, y1(i-1) + q3, y2(i-1) + r3, y3(i-1) + v3, y4(i-1) + w3);
r4 = h*f2(t(i-1) + h, y1(i-1) + q3, y2(i-1) + r3, y3(i-1) + v3, y4(i-1) + w3);
v4 = h*f3(t(i-1) + h, y1(i-1) + q3, y2(i-1) + r3, y3(i-1) + v3, y4(i-1) + w3);
w4 = h*f4(t(i-1) + h, y1(i-1) + q3, y2(i-1) + r3, y3(i-1) + v3, y4(i-1) + w3);

y1(i) = y1(i-1) +(1/6)*q1 + (1/3)*q2 + (1/3)*q3 + (1/6)*q4;
y2(i) = y2(i-1) + (1/6)*r1 + (1/3)*r2 + (1/3)*r3 + (1/6)*r4;
y3(i) = y3(i-1) + (1/6)*v1 + (1/3)*v2 + (1/3)*v3 + (1/6)*v4;
y4(i) = y4(i-1) + (1/6)*w1 + (1/3)*w2 + (1/3)*w3 + (1/6)*w4;

% q1 = h*f(t(i-1),xrg4(i-1));
% q2 = h*f(t(i-1) + h/2, xrg4(i-1) + q1/2);
% q3 = h*f(t(i-1) + h/2, xrg4(i-1) + q2/2);
% q4 = h*f(t(i-1) + h, xrg4(i-1) + q3);
%xrg4(i) = xrg4(i-1) +(1/6)*q1 + (1/3)*q2 + (1/3)*q3 + (1/6)*q4;
end
t = [a t]; % Insert t0 at the front of the t vector.
y1 = [y1a y1]; y2 = [y2a y2];y3 = [y3a y3];y4 = [y4a y4]; % Insert xa at the front of the soln vector
%% Plotting

%xexact = (K*xa)./(xa + (K-xa).*exp(-r*t));
%[T Y] = ode45(@fty, [0, 30], [4,0,0,0.2]);
%fty(2,1) = y(4);
%fty(3,1) =
%y1exact=-x(1)/(x(1)^2 + x(2)^2)^(3/2);
%y2exact = -x(2)/(x(1)^2 + x(2)^2)^(3/2);
%y1exact = -x/(x^2+y^2)^(3/2);%2*t + 2;
%y2exact = exp(t);
%fty(1,1) = y(3);
%fty(2,1) = y(4);
%fty(3,1) = -y(1)/(y(1)^2 + y(2)^2)^(3/2);
%fty(4,1) = -y(2)/(y(1)^2 + y(2)^2)^(3/2);
figure(1)
plot(t,y1,’gs’)
hold on
plot(t,y1exact,’r’)
title(‘Solution y1′)

figure(2)
plot(t,y2,’gs’)
hold on
plot(t,y2exact,’r’)
title(‘Solution y2′)

figure(3)
plot(t,y3,’gs’)
hold on
plot(t,y3exact,’r’)
title(‘Solution y3′)

figure(4)
plot(t,y4,’gs’)
%hold on
%plot(t,y1exact,’r’)
title(‘Solution y4’)

%figure(3)
%errorrg4 = xrg4 – xexact;
%semilogy(t,errorrg4)
%title(‘Error’)

%relerrorrg4 = (xrg4 – xexact)./xexact;
%figure(4)
%semilogy(t,relerrorrg4)
%title(‘Relative Error’)