# Solving $\frac{d^2y}{dx^2} + \sin(y) = 0$ with endpoint conditions of $y(0) = 0$ and $y(0.5) = 0.2618$ in MATLAB

I have this problem that i would like to solve preferably in MATLAB

Solve:

$$\frac{d^2y}{dx^2} + \sin(y) = 0$$

with endpoint conditions $y(0) = 0$ and $y(0.5) = 0.2618$.

I think i have to plot the solution. I would appreciate any kind of help.

My incomplete solutions until now is:

a=0;b=10; %limits
h=0.00001; %step

N=round((b-a)/h)+1; %nodes

x=zeros(N,1);% init solution vectro

x(1)=0; %initial condition

% y(1) = 0.2618; - i still need the second condition

for i=2:N
x(i)=x(i-1)+h*y(i-1);
y(i)=y(i-1)- h * sin(x(i-1));
end

plot(a:h:b, x);
hold on;

-

If you want a numerical solution, you can use a shooting method. It works like this. Choose $\alpha\in\mathbb{R}$ and solve by your favourite numerical method the problem $$y''+\sin(y)=0,\quad y(0)=0,\quad y'(0)=\alpha$$ on the interval $[0,1/2]$. The value found at $x=1/2$ is a function of $\alpha$ that I will denote by $s(\alpha)$. To solve the problem, you must find $\alpha$ such that $s(\alpha)=0.2618$. It is clear that $s(0)=0$. Find a value of $\alpha$ such that $s(\alpha)>0.2618$ and apply for instance the bisection method to find an approximate solution of $s(\alpha)=0.2618$.

In Mathematica you can get the solution as

u=y/.First@NDSolve[{y''[x]+Sin[y[x]]==0,y[0]==0,y[0.5]==0.2618},y,{x,0,.5}]


The command Plot[u[x],{x,0,.5},AspectRatio->Automatic] gives the plot

Finally, multiplying the equation by $y'$ and integrating once we get $$\frac12\,(y')^2+(1-\cos(y))=C,\quad C\ge0,$$ from where $$\frac{dy}{\sqrt{2\,C-(1-\cos(y))}}=\pm\,dx,$$ which can be integrated in terms of elliptic functions.

-

So i search a bit more on the internet for a matlab solution and i found this great video http://blanchard.ep.wisc.edu/bvp4cflash/bvp4cflash.html

This really explains on the same kind of problem that i have. Following his examples I tried to adjust his code to my problem and i got this (im not sure 100% is corect - but at least it plots the same plot as Julian's)

xlow = 0; xhigh = 0.5;
solinit = bvpinit(linspace(xlow,xhigh,100), [1 -1]);
sol = bvp4c(@exam1ode, @exam1bc, solinit);
xint = linspace(xlow, xhigh);
Sxint = deval(sol,xint);
plot(xint, Sxint(1,:));


exam1ode defines my ecuation

function f1 = exam1ode(x, y)
f1 = [y(2) -sin(y(1))];
end


exam1bc defines the boundary conditions

function f2 = exam1bc(ya, yb)
f2 = [ya(1) yb(1)-0.2618];
end

-