Midpoint numerical integration without a built in function

I need some help building a matlab script to solve dy/dt = y*t^3-1.5*y using the midpoint method. I have solved this using Euler's and the below code

a = 0;
b = 2;

h=.5;
T = 0:h:2;

y = zeros(1,((b-a)/h+1));
y(1) = 1;
phi1 = y.*T.^3-1.5.*y;
for i = 2:length(T)
    phi1 = y(i-1).*T.^3-1.5.*y(i-1);
    y(i) = y(i-1) + phi1(i-1).*h;
end
plot(T,y,':bo')

But solving cannot figure out the midpt method as I know the +1/2 intervals are tough on MATLAB.

Below is what I have for midpoint and I know it is very wrong.

thanks in advance.

a = 0;
b = 2;
h= .5;
h2 = h/2;
t = 0:h2:2;

y = zeros(1,9); 
y(2) = 1;
phi3 = y.*t.^3-1.5.*y;
for i = 3:(length(t))
    Y2 = y(i-2)+ (y(i-2).*(t).^3-1.5.*y(i-2))*(h/2);
    phi3 = Y2.*t.^3-1.5.*Y2;
    y(i) = y(i-2) + phi3(i-1).*h;
end

plot(t,y);

[Roger Stafford]

It pains me to see you use numerical integration for this differential equation when its solution can be expressed so easily as

   y = K*exp(t^4/4-1.5*t)

where K depends on initial conditions.

Answer 1

[Geoff Hayes]

Aggie - the midpoint method should be very similar to your Euler implementation, with just a couple of minor changes (for example the step size). The above code gets a little confusing with the re-calculation of phi3 at every iteration (as it is vector of which you only use one element from), so you may want to try an alternative approach where you define a function handle to your equation and evaluate that instead

 F = @(t,y)y*t^3-1.5*y;

Then, the difference relation for the midpoint algorithm can be written as

 y(n) = y(n-1) + h*F(t(n-1) + h/2, y(n-1) + (h/2)*F(t(n-1),y(n-1)));

So your code becomes

 a = 0;
 b = 2;
 h = 0.5;
 F = @(t,y)y*t^3-1.5*y;

 t = a:h/2:b;
 y = zeros(size(t));

 % set the initial condition
 y(1) = 1;

 for n=2:length(t)
      y(n) = y(n-1) + h*F(t(n-1) + h/2, y(n-1) + (h/2)*F(t(n-1),y(n-1)));
 end

 plot(t,y,':ro');

Note that the body of the for loop has been reduced to one line, and is quite different from yours. Is there a reason that you started iterating at i equal to 3, and referred to i-2 in your equation?

[Konner Brickey]

Hello,

There is a slight problem with this code. The line

t = a:h/2:b;

should be

t = a:h:b