// Режим отображения
mode(0);

// Отображать предупреждение для исключений при выполнении операций с плавающей запятой
ieee(1);

a1 = 10;
a2 = 10;
L = 10;
dx = 0.05;

t1 = sqrt(L/((1+a1/a2)*a1));
t2 = (a1*t1)/a2;
Xm = (a1*t1)*t1;

X = 0:dx:L;
n = max(size(X));
T = zeros(1,n);
Y = zeros(1,n);

for i = 2:1:n
  if X(i)<=Xm then
    a = a1/2;
    b = Y(i-1);
    c = -dx;
    D = b*b-(4*a)*c;
    dt1 = (-b+sqrt(D))/(2*a);
    dt2 = (-b-sqrt(D))/(2*a);
    if dt2>0 then
      if dt1<dt2 then
        dt = dt1;
      else
        dt = dt2;
      end;
    else
      dt = dt1;
    end;
    T = mtlb_i(T,i,T(i-1)+dt);
    Y = mtlb_i(Y,i,Y(i-1)+a1*dt);
  else
    a = -a2/2;
    b = Y(i-1);
    c = -dx;
    D = b*b-(4*a)*c;
    dt1 = (-b+sqrt(D))/(2*a);
    dt2 = (-b-sqrt(D))/(2*a)
    if dt2>0 then
      if dt1<dt2 then
        dt = dt1;
      else
        dt = dt2;
      end;
    else
      dt = dt1;
    end;
    T = mtlb_i(T,i,T(i-1)+dt);
    Y = mtlb_i(Y,i,Y(i-1)-a2*dt);
  end;
end;

X1 = X;
Y1 = Y;
leg1 = "t = "+msprintf("% 10.4f",t1+t2)+"s; a1 = "+msprintf("% 10.4f",a1)+"; a2 = "+msprintf("% 10.4f",a2)

a1 = 1;
a2 = 3;
L = 100;
dx = 0.05;

t1 = sqrt(L/((1+a1/a2)*a1));
t2 = (a1*t1)/a2;
Xm = (a1*t1)*t1;

X = 0:dx:L;
n = max(size(X));
T = zeros(1,n);
Y = zeros(1,n);

for i = 2:1:n
  if X(i)<=Xm then
    a = a1/2;
    b = Y(i-1);
    c = -dx;
    D = b*b-(4*a)*c;
    dt1 = (-b+sqrt(D))/(2*a);
    dt2 = (-b-sqrt(D))/(2*a);
    if dt2>0 then
      if dt1<dt2 then
        dt = dt1;
      else
        dt = dt2;
      end;
    else
      dt = dt1;
    end;
    T = mtlb_i(T,i,T(i-1)+dt);
    Y = mtlb_i(Y,i,Y(i-1)+a1*dt);
  else
    a = -a2/2;
    b = Y(i-1);
    c = -dx;
    D = b*b-(4*a)*c
    dt1 = (-b+sqrt(D))/(2*a)
    dt2 = (-b-sqrt(D))/(2*a)
    if dt2>0 then
      if dt1<dt2 then
        dt = dt1
      else
        dt = dt2
      end;
    else
      dt = dt1;
    end;
    T = mtlb_i(T,i,T(i-1)+dt);
    Y = mtlb_i(Y,i,Y(i-1)-a2*dt);
  end;
end;

X2 = X;
Y2 = Y;
leg2 = "t = "+msprintf("% 10.4f",t1+t2)+"s; a1 = "+msprintf("% 10.4f",a1)+"; a2 = "+msprintf("% 10.4f",a2)

a1 = 4.9;
a2 = 4.9;
L = 100;
dx = 0.05;

t1 = sqrt((2*L)/((1+a1/a2)*a1));
t2 = (a1*t1)/a2;
Xm = (a1*t1)*t1;

X = 0:dx:L;
n = max(size(X));
T = zeros(1,n);
Y = zeros(1,n);

for i = 2:1:n
  if X(i)<=Xm then
    a = a1/2;
    b = Y(i-1);
    c = -dx;
    D = b*b-(4*a)*c;
    dt1 = (-b+sqrt(D))/(2*a)
    dt2 = (-b-sqrt(D))/(2*a)
    if dt2>0 then
      if dt1<dt2 then
        dt = dt1;
      else
        dt = dt2;
      end;
    else
      dt = dt1;
    end;
    T = mtlb_i(T,i,T(i-1)+dt);
    Y = mtlb_i(Y,i,Y(i-1)+a1*dt);
  else
    a = -a2/2;
    b = Y(i-1);
    c = -dx;
    D = b*b-(4*a)*c
    dt1 = (-b+sqrt(D))/(2*a)
    dt2 = (-b-sqrt(D))/(2*a)
    if dt2>0 then
      if dt1<dt2 then
        dt = dt1;
      else
        dt = dt2;
      end;
    else
      dt = dt1;
    end;
    T = mtlb_i(T,i,T(i-1)+dt);
    Y = mtlb_i(Y,i,Y(i-1)-a2*dt);
  end;
end;

X3 = X;
Y3 = Y;
//leg3 = ['t = ', num2str(t1+t2,'% 10.4f'), 's; a1 = ', num2str(a1,'% 10.4f'), '; a2 = ',num2str(a2,'% 10.4f')]

plot(X1,Y1,X2,Y2,X3,Y3);
//grid on;
//xlabel('x(t)');
//ylabel('y(t)');

//legend(leg1,leg2,leg3);