function dummy = soliton(n, ct, flow)

format long e;
deltax = 0.001;
xmax = 5;
x = 0;
y = 0;
yprime = 1;
ypprime = 0;


while x(length(x)) <= xmax & y < pi/2
    
    
    
    
    xi = x(length(x)) + deltax;
    x = [x, xi];
    yprimei = yprime(length(yprime))+ deltax*ypprime(length(ypprime));
    yprime =[yprime, yprimei];
    yi = y(length(y)) + deltax*yprimei;
  
    
    y = [y, yi];
    switch flow
        case 'Ricci'
            ypprimei = (-ct*xi*yprimei)/(cos(yi)^2) - ((n-2)*yprimei)/(xi*cos(yi)^2) + ((n-2)*tan(yi))/(xi^2) - (yprimei^2/tan(yi)) + yprimei/xi;
            
            
        case 'Mean'
            ypprimei = (-ct*xi*yprimei)/(cos(yi)^2) - ((n-1)*yprimei)/(xi*cos(yi)^2) + ((n-1)*tan(yi))/(xi^2);
    
        case 'Gauss'
            ypprimei = (-ct*(xi)^n*yprimei)/((sin(yi)^(n-1))*(cos(yi)^2)) + ((n-1)*yprimei)/(xi) - ((n-1)*yprimei^2)/(tan(yi));
    
        case 'Yamabe'
            ypprimei = (-ct*xi*yprimei)/((n-1)*cos(yi)^2) - ((n-2)*yprimei)/(2*(cos(yi)^2)*xi) + ((n-2)*tan(yi))/(xi^2) - (yprimei^2)/tan(yi) - (yprimei*(n-4))/(2*xi);
    
    
        case 'Einstein'
            ypprimei = (-ct*xi*yprimei)/((n-2)*cos(yi)^2)- (yprimei*(n-3))/(2*xi)+ ((n-2)*tan(yi))/(xi^2)- (yprimei^2)/tan(yi)-((n-3)*yprimei)/(2*(cos(yi)^2)*xi);
        otherwise
            disp('Unrecognized Flow, options are Ricci, Mean, Gauss, Yamabe, or Einstein')
            break
    end
    ypprime = [ypprime, ypprimei];

    
    
end

bprime = 0;
b = 0;
for i = 1:length(x)
    bprime(i) = tan(y(i));
    b(i) = trapz(bprime)*deltax;
end
plot (0,0)
hold on
%plot(x, bprime)
%plot(-x, -bprime)
plot(b, x)
plot(b,-x)


hold off
axis ([0 2*xmax -xmax xmax])
%axis([-xmax xmax -xmax xmax])
axis square;

end
