Solving DE using taylor's method

%Program to solve ordinary differential equations using modified Taylor’s method
clc;
clear;
close all;
syms x y;
%Define function
f(x,y)=1/(x^2+y);
y0=4;
h=0.1;
x0=4;
xn=4.1;
% Define the order of the Taylor series expansion
n=1;
y(1)=y0;
y(2)=f(x,y);
% Compute the values of h^n/n! for n = 0, 1, 2, 3,...,
ht=h.^(0:n)./factorial(0:n);
for i=2:n
% To compute the values of derivatives
y(i+1)=diff(y(i),x)
end
for i=1:n
x=x0+i*h;
c=eval(y);
% Compute the values of yb using the Taylor series expansion
yb=sum(c.*ht);
fprintf('The value of y(%f) is %f\n',x,yb);
end