Solving a problem with the Shooting Newton Method in Matlab

I tried to implement the Shooting method using Newton method to solve a differential equation, but the exact solution recommended by the professor does not coincide with the numerical solution.

Equation:

-y''+exp(-x)/2(y'^2+y^2)=exp(-x)/2*cos(x)+2sin(x)

BC:

2y(0)-y'(0)=0

y(pi)=exp(pi)

Exact solution:

y=exp(x)+sin(x)

The graphs should coincide, but they don't.

Can you help me?

For testing

>> [x,u]= NonLinearShootingNewtonExercise(20,1,0.0001);

>> y=exp(x)+sin(x);

>> plot(x,u,'r*',x,y,'k-')

My code is:

[x,u]= NonLinearShootingNewtonExercise(20,1,0.0001);
y=exp(x)+sin(x);
plot(x,u,'r*',x,y,'k-')

function [x,u] = NonLinearShootingNewtonExercise(N,sk,tol)  
%Shooting method using Netwod method

%Equation:
%-y''+exp(-x)/2(y'^2+y^2)=exp(-x)/2*cos(x)+2sin(x)
%BC:
%2y(0)-y'(0)=0
%y(pi)=exp(pi)

%Exact solution:
%y=exp(x)+sin(x)

%Input
%N: number of discretization point
%s0: initial value of s
%tol: tollerance

% Output:
% x: Vector containing the points at which the solution is approximated
% u: Approximated solution

%For testing
% [x,u]= NonLinearShootingNewtonExercise(20,1,0.0001);
% y=exp(x)+sin(x);
% plot(x,u,'r*',x,y,'k-')

epi=exp(pi);
h=pi/N;
x=0:h:pi;
ds=1+tol;

while ds>=tol
    [xx,uvUV]=ode23(@fun,x,[sk 2*sk 1 2]); %I have 4 initial values (u,v,U,V)
   %uvUV is an array of 4 coloumns and N+1 rows
   
   phik=uvUV(N+1,1)+uvUV(N+1,2)-epi;
   dphik=uvUV(N+1,3)+uvUV(N+1,4); %derivate of phi
   
   skp1=sk-phik/dphik;
   ds=abs(skp1-sk);
   sk=skp1;
   
end

u=uvUV(:,1);
end

function uvp= fun(x,uvUV)
uvp=[uvUV(2); (exp(-x)/2*(uvUV(2)^2+uvUV(1)^2))-(exp(-x)/2+cos(x)+2*sin(x));...
    uvUV(4);(exp(-x)*uvUV(1)*uvUV(3))+(exp(-x)*uvUV(2)*uvUV(4))];
%u is the first component and v is the second component
end

[Francesca]

%Esercizio n.1 del 22 ottobre

function [x, u] = NonLinearShootingNewtonExercise(N, s0, tol)  
% Shooting method using Newton's method
% Equation:
% -y''+exp(-x)/2(y'^2+y^2)=exp(-x)/2+cos(x)+2sin(x)
% BC:
% 2y(0)-y'(0)=0
% y(pi)=exp(pi)
% Exact solution:
% y=exp(x)+sin(x)
% Input
% N: number of discretization points
% s0: initial value of s
% tol: tolerance
% Output:
% x: Vector containing the points at which the solution is approximated
% u: Approximated solution

% For testing
% [x, u] = NonLinearShootingNewtonExercise(20, 1, 0.0001);
% y = exp(x) + sin(x);
% plot(x, u, 'r*', x, y, 'k-')

esponenzialepi = exp(pi);
h = pi / N;
x = 0:h:pi;
ds = 1 + tol;
sk = s0;
while ds >= tol
    [xx, uvUV] = ode23(@fun, x, [sk 2*sk 1 2]);
    phik = uvUV(N+1, 1)  - esponenzialepi;
    dphik = uvUV(N+1, 3) ; % derivative of phi
    skp1 = sk - phik / dphik;
    ds = abs(skp1 - sk);
    sk = skp1;
end
u = uvUV(:, 1);
u(end) = exp(pi); % Set the boundary condition at x=pi
end

function uvp = fun(x, uvUV)
u = uvUV(1);
v = uvUV(2);
du = uvUV(3);
dv = uvUV(4);

uvp = [v; ...
    (exp(-x) / 2 * (v^2 + u^2)) - exp(-x) / 2 - cos(x) - 2 * sin(x); ...
    dv; ...
    (exp(-x) * u * du) + (exp(-x) * v * dv)];
end
 

Answer 1

[Torsten]

sol = fsolve(@fun,23,optimset('TolFun',1e-12,'TolX',1e-12));

Equation solved, inaccuracy possible.

The vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

fun(sol)

ans = -1.1977e-12

[X,Y] = ode_solver(sol);
figure(1)
hold on
plot(X,Y(:,1))
plot(X,exp(X)+sin(X))
hold off
grid on

sol = fsolve(@fun,30,optimset('TolFun',1e-12,'TolX',1e-12));

Equation solved, solver stalled.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared and the vector of function values
is near zero as measured by the value of the function tolerance.

fun(sol)

ans = 3.0909e-13

[X,Y] = ode_solver(sol);
figure(2)
hold on
plot(X,Y(:,1))
plot(X,exp(X)+sin(X))
hold off
grid on

function res = fun(p)
  [X,Y] = ode_solver(p);
  res = 2*Y(end,1)-Y(end,2);
end

function [X,Y] = ode_solver(p);
  fun_ode = @(x,y)[y(2);-exp(-x)/2*cos(x)-2*sin(x)+exp(-x)/2*(y(2)^2+y(1)^2)];
  y0 = [exp(pi);p];
  tspan = [pi,0];
  [X,Y] = ode15s(fun_ode,tspan,y0,odeset('RelTol',1e-12,'AbsTol',1e-12));
end

[Francesca]

How can I modify the code while maintaining the structure and without changing my commands too much?