Question

Write Matlab code that solves (dy/dt)= yt^3 - 1.5y on the interval t= 0 to 2 where y(0) = 1. As well as (dy/dx) = (1+4x)y^.5 on the interval x = 0 to 1 where y(0) = 1 using Matlab functions that solve any ODE using the below method

• Runge-Kutta Method 4th Order

Functions need to call:

• Error desired
• number of iterations
• Step size
• Initial Interval
• y(0)

And return

• Result
• Error of that result

Answer 1

`Hey,

Note: In case of any queries, just comment in box I would be very happy to assist all your queries

All previous parts as well as current parts combined

clc
clear all
close all
f=@(t,y) y*t^3-1.5*y;
g=@(x,y) (1+4*x)*y^(0.5);
disp('Equation 1 using ode45');
[T,Y]=ode45(f,[0,2],1)
plot(T,Y);
title('Equation 1 plot ode45');
disp('Equation 2 using ode45');
[X,Y]=ode45(g,[0,1],1)
figure;
plot(X,Y);
title('Equation 2 plot ode45');
figure;
disp('Equation 1 using euler');
[T,Y]=eulerSystem(f,[0,2],1,0.05)
plot(T,Y);
title('Equation 1 plot euler');
disp('Equation 2 using euler');
[X,Y]=eulerSystem(g,[0,1],1,0.05)
figure;
plot(X,Y);
title('Equation 2 plot euler');
figure;
disp('Equation 1 using heun');
[T,Y]=heun(f,[0,2],1,0.05)
plot(T,Y);
title('Equation 1 plot heun');
disp('Equation 2 using heun');
[X,Y]=heun(g,[0,2],1,0.05)
figure;
plot(X,Y);
title('Equation 2 plot heun');
figure;
disp('Equation 1 using midpoint');
[T,Y]=midpoint(f,[0,2],1,0.05)
plot(T,Y);
title('Equation 1 plot midpoint');
disp('Equation 2 using midpoint');
[X,Y]=midpoint(g,[0,1],1,0.05)
figure;
plot(X,Y);
title('Equation 2 plot midpoint');

figure;
disp('Equation 1 using runge4');
[T,Y]=runge4(f,[0,2],1,0.05)
plot(T,Y);
title('Equation 1 plot runge4');
disp('Equation 2 using runge4');
[X,Y]=runge4(g,[0,1],1,0.05)
figure;
plot(X,Y);
title('Equation 2 plot runge4');
function [t,y] = midpoint(f,tspan,ya,h)
a=tspan(1);
b=tspan(2);
n = floor((b - a) / h);
halfh = h / 2;
y(1,:) = ya;
t(1) = a;
for i = 1 : n
t(i+1) = t(i) + h;
z = y(i,:) + halfh * f(t(i),y(i,:));
y(i+1,:) = y(i,:) + h * f(t(i)+halfh,z);
end
end
function [t,y] = heun(f,tspan,ya,h)
a=tspan(1);
b=tspan(2);
n = floor((b - a) / h);
halfh = h / 2;
y(1,:) = ya;
t(1) = a;
for i = 1 : n
t(i+1) = t(i) + h;
g = f(t(i),y(i,:));
z = y(i,:) + h * g;
y(i+1,:) = y(i,:) + halfh * ( g + f(t(i+1),z) );
end
end

function [x,y]=runge4(f,tspan,y0,h)
x = tspan(1):h:tspan(2); % Calculates upto y(3)
y = zeros(length(x),1);
y(1,:) = y0; % initial condition
for i=1:(length(x)-1) % calculation loop
k_1 = f(x(i),y(i,:));
k_2 = f(x(i)+0.5*h,y(i,:)+0.5*h*k_1);
k_3 = f((x(i)+0.5*h),(y(i,:)+0.5*h*k_2));
k_4 = f((x(i)+h),(y(i,:)+k_3*h));
y(i+1,:) = y(i,:) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h; % main equation
end
end

function [t,y]=eulerSystem(Func,Tspan,Y0,h)
t0=Tspan(1);
tf=Tspan(2);
N=(tf-t0)/h;
y=zeros(length(Y0),N+1);
y(:,1)=Y0;
t=t0:h:tf;
for i=1:N
y(:,i+1)=y(:,i)+h*Func(t(i),y(:,i));
end
end

Kindly revert for any queries

Thanks.