Runge-Kutta to solve differential equations.

%Runge-Kutta to solve differential equations.
%Author: minster
% solve x’ = f(t,x); x(a) = xa
clear all
close all

%%Initialize
a = 0; %initial time
b = 10; %final time
xa = 1; %initial solution
%r = 0.5; %parameters for the DE below
%K = 10;
f = inline(‘0.5*(1-x/10)*x’,’t’,’x’);
r = 0.5; %parameters for the DE belowk
k = 10;
n=100; %Number of time steps
h = (b-a)/n; %step size
t = (a+h:h:b); %vector of times
xrg4 = t; %solution vector. Copy time into it to create it with appropriate length.

%% Begin Runge-Kutta 2nd order and 4th order.
F1 = h*f(a,xa); %
F2 = h*f(a + h/2, xa + F1/2);
F3 = h*f(a + h/2, xa + F2/2);
F4 = h*f(a + h, xa + F3);
xrg4(1) = xa +(1/6)*F1 + (1/3)*F2 + (1/3)*F3 + (1/6)*F4;

for i = 2:n
tk = a + k*h;
F1 = f(tk, yk);
F2 = f(tk+(h/2), yk + (h/2)*F1);
F3 = f(tk+h/2, yk + (h/2)*F2);
F4 = f(tk+h, yk+h*F3);
yk = yk + (h/6)*(F1 + 2*F2 + 2*F3 + F4);
y(:,k+1) = yk;

% 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)*F1 + (1/3)*F2 + (1/3)*F3 + (1/6)*F4;
end
t = [a t]; % Insert t0 at the front of the t vector.
xrg4 = [xa xrg4]; % Insert xa at the front of the solution vector, x.

%% Plotting

xexact =(K*xa)./(xa + (K-xa).*exp(-r*t));

figure(2)
plot(t,xrg4,’gs’)
hold on
plot(t,xexact,’r’)
title(‘Solution Curves’)

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

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