clear;
R=2.8e-1;L=1e-5;E=126;T=1/12000;Alpha=9.9999;
Theta1=R/L;
Theta2=E/L;
Phi1=exp(-Theta1*T);
Phi2=Theta2/Theta1;
Phi3=Theta1*T;
Xmax=E/R;
for i=1:1:999,
u(i)=i/500-0.999;
u1=u(i);
ter1=1+u1-u1*sign(u1);
ter2=1-u1*sign(u1);
  if ter1==0,
     x(i)=Xmax;
  else
     x(i)=Xmax*ter1/ter2;
  end
end
plot(u,x,'r-'),axis([-1 1 -Xmax*3 Xmax*5]);
pause;
%
xr(1)=x(1);
for i=1:1:998,
	Mua=u(i);
	Ter1=(Phi1^(1-Mua)-Alpha);
	Ter2=(Phi2/2)*(Phi3*Phi1^(1-Mua)*sign(Mua)*Mua+2);
        Ter3=(Phi2/2)*Phi3*sign(Mua)*u(i+1);
       xr(i+1)=-(Ter2+Ter3+Alpha*x(i)-x(i+1))/Ter1;
end
clg;
plot(u,x,'r-',u,xr,'g-'),axis([-1 1 -Xmax*3 Xmax*5])
pause;
ur(1)=u(1);
for i=1:1:998,
	Mua=ur(i);
	Ter1=(Phi1^(1-Mua)-Alpha)*xr(i);
	Ter2=(Phi2/2)*(Phi3*Phi1^(1-Mua)*sign(Mua)*Mua+2);
        Ter3=(Phi2/2)*Phi3*sign(Mua);
       ur(i+1)=-(Ter1+Ter2+Alpha*x(i)-x(i+1))/Ter3;
end
clg;
plot(u,x,'r-',u,xr,'g-',ur,xr,'b-'),axis([-1 1 -Xmax*3 Xmax*5])

 