dx/dt=y,dy/dt=-sinx,求大神帮忙编一个MATLAB的程序，用龙格库塔法解这个方程组，求关于x,y的数值解
建立m文件：
function dx=dfun(t,x) %函数名为dfun，参数为t与x
dx=[x(2);-sin(x(1))]; %以向量形式表示方程
输入：
clear
ts=-15:0.05:15; %步长取0.05
x0=[1,0]; %设定参数初值
options=odeset('reltol',1e-6,'abstol',1e-9); %提高精度
[t,x]=ode45(@dfun,ts,x0,options); %调用ode45计算
plot(x(:,1),x(:,2)),grid %作出y(x)图形
axis equal
gtext('\fontsize{12}x'),gtext('\fontsize{12}y') %标记字体x
但以上并非曲线y=f(x)的完整形状（调整ts的范围也无济于事），原因是y为x的周期函数，而数值解只能求出初值附近的解
本题可以求出y=f(x)的解析表达式
由dx/dt=y，dy/dt=-sinx，得
dy/dx=(dy/dt)*1/(dx/dt)=-sinx/y
分离变量，积分得
y^2=2*cos(x)+C，其中C为常数
代入初始条件y(1)=0，可求得C=-2*cos(1)
∴y^2=2*cos(x)-2*cos(1)，此式为原方程组的解析解
利用ezplot命令可绘制出完整图像
clear
syms x y
ezplot(y^2-2*cos(x)+2*cos(1),[-8,8,-3,3])
axis equal
axis([-8,8,-3,3])
grid on
另外，改变初值将得到不同的图形（为什么？请思考），例如
初值改为：x=1,y=√[2*(cos(1)+1)]-10^(-5)
初值改为：x=1,y=√[2*(cos(1)+1)]+10^(-5)
