

clear all

PSI=linspace(0.1,1,1001);

tb=8/9*sqrt(3);

T=[0,1/2,1,2,3]*tb;

for ia = 1:length(T);
    
    t=T(ia);
    
    X1=t*PSI+sqrt(1./PSI-1);
    X2=t*PSI-sqrt(1./PSI-1);
    
    figure(100)
    if ia==5;
        ip=1;
    else;
        ip=ia
    end;
    subplot(2,2,ip)
    plot(X1,PSI,'r',X2,PSI,'b','linewidth',2);
   
    
end;


XI=linspace(-1,4,101);

F=1./(1+XI.^2);

T=linspace(0,3*tb,101);

for ia=1:length(XI)
    
    f=F(ia);
    X=f*T+XI(ia);
    figure(500)
    plot(X,T,'b','linewidth',0.5)
    hold on
end;


xip=1/sqrt(3)*linspace(0.1,7,501);


tp=((1+xip.^2).^2)./(2*xip);
xp=(1./(1+xip.^2)).*tp+xip;
    figure(500)
    plot(xp,tp,'g','linewidth',3)
    hold off

    
    
    
    
    
    





