Hello, I provide the code for explicit and implict theta method

I tried to do it but I got an error.

please apply theta method for theta = , / , with = . for y'=x-y^2 on [0,1] with 0.01 step size and then do the following.

a) Plot in one figure the numerical solutions together with the exact solution.

b) Plot in one figure the absolute global errors || = |( ) | for = 0, 1/2, 1.

c) attached the code after you modify it

thanks

function [T , x_node, y_num, y_exact, G_Error] = ThetaMethodDE(a,b,h,theta,M,f,y,y0,tolerance,option)

% % Theta methods for solving ODE using iteration method to solve the non-linear equation as

% % sub-problem which is defined in Numerical Solution of Ordinary Differential Equations

% % By Kendall E. Atkinson, Weimin Han, and David E. Stewart 2009 John Wiley & Sons, Inc.

% % ISBN 978-0470-04294-6

% % Inputs:

% % a: the intial point in interval [a,b]

% % b: the terminal point in interval [a,b]

% % h: stepsize of subinterval

% % theta: value of a parameter theta

% % tolerance: the stopping criteria for iteration method

% % M : number of different h

% % f: the gradiant of y as y' = f(x1,x2)

% % y0: initial condition at a x0

% % y: exact solution

% % option : 'true' display table of outputs

% % 'false' hide table of outputs

% % Outputs: T: table of numerical solutions at the points x_i, i=0,...,n ,

% % x_node: Array of mesh points

% % y_num: Array of numerical values at x_node

% % y_exact: Array of exact solutions at x_node

% % G_Error: Array of absolute globel error e_n = |(y(x_n)-y_n)|

%% the exact solution

z=2/3 * x.^(3/2);

rootx=x.^(1/2);

sqrt3=sqrt(3);

bk1=sqrt3 *besselk(-2/3 , z)/pi;

bi1= besseli(-2/3 , z);

bk2= sqrt3 *besselk(1/3 , z)/pi;

bi2=besseli(1/3, z);

ye=-rootx .* (bk1 - bi1) ./(bk2+bi2);

%% the y derivative = x-y^2 on the interval [0,1] & intial condition y(o)=0

StepSize = zeros(1,M);

for k = 1:M

n= (b-a)/h;

x_node = linspace(a,b,n+1);

y_num = zeros(1,n+1);

y_num(1) = y0;

for i = 1:n

if theta == 0

y_num(i+1) = y_num(i) + h*f(x_node(i),y_num(i));

else

fun_val = f(x_node(i),y_num(i));

expEuler = y_num(i) + h*fun_val;

count = 0;

Iter_diff = 1;

while Iter_diff > tolerance && count < 10

y_num(i+1) = y_num(i) + h*((1-theta)*fun_val + theta*f(x_node(i+1),expEuler));

Iter_diff = abs(expEuler-y_num(i+1));

expEuler = y_num(i+1);

count = count +1;

end

end

end

% computing absolute global error e_n =|y(x_n)- y_n|

y_exact = y(x_node);

G_Error = abs(y_exact - y_num);

% Table of Outputs

switch(option)

case('true')

Ttitle = string( [' ' ...

' Output for theta = ' num2str( theta ) ] );

disp(Ttitle)

T = unique(array2table([x_node',y_num',y(x_node)', G_Error'], ...

'VariableNames', ...

{'x nodes','Numerical Value','Exact','|Error|'}));

display(T)

case('false')

T =[];

end

%Update stepsize

StepSize(k) = h;

h = h/2;

end