% Calculate dy/dx = F(y,x), y(a) = y0, a <= x <= b
% using  Midpoint method

clc;                                               % Clears the screen
clear all;
a = 0;          % a is initial value of x
b= 1;           % b is final value of x
y0=2;             % y(a) = y0;
n=10;           % n = number of time steps
h = (b-a)/n;    % h = time step size
x = [a:h:b];    % x = [a, a+h, a+2h, ..., b]                                          % step size
y = zeros(size(x));
y(1) = y0;      % initial condition set

F_xy = @(t,r) t*r;    % F_xy is a function of two variables
                    % the first variable is the independent variable and
                    % the second is the dependent variable of the ode
                    
for i=1:n                              % calculation loop
    k_1 =  F_xy(x(i),y(i));
    k_2 =  y(i) +  h*k_1;
    y(i+1) = y(i) + 0.5*h*(k_1 + F_xy(x(i+1),k_2));  % main equation
end

figure(3);
plot(x,y,'o', x, y0*exp(0.5*(x-a).^2),'+'); grid
xlabel('x')
ylabel('y')
title('Trapezoid solution to dy/dx = F(x,y)');   % change this as necessary for 
                            % different right hand sides