clear
to = 1;
t = 160;
t1=t;
L=0.5;
N = 60;


x(1:2*N) = 0;
options = odeset('RelTol',1e-7,'OutputFcn','odeplot','OutputSel',[3 6]);
[t x] = ode45('fputest1', [0:L:t], x, options, N);


for H=to/L:t1/L, EE(H)=0; hh(H)=H*L; hhh=H*L; 
        for J=1:N-1, 
        YY(H,J)=x(H,J);
        EE(H)=EE(H)+x(H,N+J)^2/2+(x(H,J+1)-x(H,J))^2+(x(H,J+1)-x(H,J))^4/2;
end
end