在床上躺了一会儿，很快就明白过来我的随机数不太对了……

reproduct
y=100

修正了一下初速度方向

细节

来十万个

哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈

clc
clear all
f=dsolve('cos(t)*(sin(pi*t/T)^2)==D2y');
syms T theta R y0 y0d x0 x0d t x y
s=diff(f,t);
f0=subs(f,T,35);
syms C5 C4
C5=solve(subs(f0,t,0)==y0,C5);
df=diff(f,t);
dff=subs(df,T,35);
C4=solve(subs(dff,t,0)==y0d,C4);
f1=eval(f);
f2=x0+t*x0d;
f3=subs(f1,t,(x-x0)/x0d);
s0=subs(eval(s),t,(x-x0)/x0d);
u=subs(s0,T,35);
u=subs(u,y0,R*sin(theta));
u=subs(u,x0,R*cos(theta));
u=subs(u,R,1);
f4=subs(f3,y0,R*sin(theta));
f4=subs(f4,x0,R*cos(theta));
f4=subs(f4,T,35);
f4=subs(f4,R,1);
d=1000;
c=0;
h=1;
for i=1:100
x0d=2*(rand-0.5);
y0d=2*(rand-0.5);
theta=rand*pi;
g=eval(f4);
u0=eval(u);
if atan(abs(cot(theta)))+atan(abs(y0d/x0d))<=pi/2
c(h)=solve(g==d);
m(h)=subs(eval(u),x,c(h));
XD(h)=x0d;
h=h+1;
end
end
cmin=min(c);
cmax=max(c);
k=linspace(-10000,10000,1000);
cc=hist(c,k);
cc=cc/length(c);
bar(k,cc)
dd=hist3(XD,m,k);
bar3(XD,m)

http://www.nicovideo.jp/watch/sm17727986

这使得通过使用组合了相位同步延迟的孤立阿秒紫外脉冲或脉冲序列的强烈的少量光学周期的红外激光脉冲分别测量作为在红外光和紫外脉冲或脉冲序列之间的可调节延迟τ的函数的光电子产量的条纹光电子光谱和双光子双路径光电发射干涉图像的记录成为可能

㔔

function [Y] = RK45(t,X,f,h)
K1=f(t,X);
K2=f(t+h/2,X+h/2*K1);
K3=f(t+h/2,X+h/2*K2);
K4=f(t+h,X+h*K3);
Y=X+h/6*(K1+2*K2+2*K3+K4);
end
以上是4阶龙格库塔法的代码:
自己写函数，存为f.m
function dxdt = f (t,x)
dxdt(1)=exp(x(1)*sin(t))+x(2);
dxdt(2)=exp(x(2)*cos(t))+x(1); % x(1)是你的f，x(2)是你的g
dxdt=dxdt(:);
end
自己给出t0,x0,h的值（初始时间，初值，步长）
如果求t0到t1的轨迹的话:给个例子如下
t0=0;t1=5;h=0.02;x0=[-1;-1];
T=t0:h:t1;X=zeros(length(x0),length(T));X(:,1)=x0;
for j=1:length(T)-1
X(:,j+1)=RK45(T(j),X(:,j),@(t,x) f(t,x),h);
end
plot(T,X(1,:));
hold on;
plot(T,X(2,:),'r');
具体参数自己设置

https://pan.baidu.com/s/1aFEMBTvhQRjdkZtIRzE8uw