%% Solution Assignment 1 (MATH246, Spring 2024)
%% Problem 1(a): y'=\frac{y}t+1,\;y(1)=2
% This is a linear ODE:
%     y'-\frac1t y=1 with a(t)=-\frac1t,\;A(t)=-\ln t (assuming t>0), multiply by integrating factor e^{A(t)}=\frac1t:
%     \frac{d}{dt}\left[\frac1t y(t)
% \right]
% = \frac1t integrate:
%     \frac1t y(t) = \ln t + C (assuming t>0)
% initial condition t_0=1, y_0=2 gives 2 = C and
%     y(t) = t\cdot \ln t+2t
% This solution exists for t\in (0,\infty). This justifies in retrospect the above assumptions.
%% Problem 1(b): y'=\frac{y^2}t and (i) y(1)=2, (ii) y(1)=0
% This is a separable ODE y'=f(t)g(y) with g(y)=y^2.
% Find all y such that g(y)=0: y_*=0 is stationary (constant) solution.
% Separate variables  \frac1{g(y)}dy=f(t)dt and integrate:
%     \int y^{-2}dy=\int t^{-1}dt gives -y^{-1} = \ln t+C (assuming t>0).
% initial condtion(i): t_0=1 and y_0=2 gives -\frac12 = C, hence
%      y(t)=\frac{1}{\frac12-\ln t}
% This solution exists for t\in (0,\sqrt e).
% initial condition (ii): we get the constant solution y(t)=0, solution exists for t\in (0,\infty)
% interval (-\infty,\infty) is also okay, but it is not clear what ODE means for t=0.
%% Problem 1(c): y'=-\frac{2t+y}{t+2y}, \; y(1)=2
% We can write this as (2t+y)dt+(t+2y)dy = 0, hence M(t,y)=2t+y,\;N(t,y)=t+2y.
% Check whether this is exact ODE: we need \partial_y M = \partial_t N
% We have \partial_yM = 1 and \partial_t N=1, hence ODE is exact, and we can find H(t,y) with \partial_t H = M, \; \partial_y H=N.
% \partial_t H(t,y) = M = 2t+y gives H(t,y)=t^2+ty+g(y) where we still have to find the function g(y).
% \partial_yH(t,y)=N=t+2y gives t+g'(y)=t+2y \iff g'(y)=2y (note that this does not depend on t), hence g(y)=y^2 (no need for constant, we just need to find one function H). We obtain
%      H(t,y)=t^2+ty+y^2
% Note: we could also first use \partial_yH=N and then \partial_tH=M, this gives the same result. Just "guessing" H(t,y) and then showing that it works is not acceptable.
% We obtain the solution in implicit form t^2+ty+y^2=C. The initial condition gives t_0=1,\;y_0=2 gives 1+2+4=C.
% Solution of IVP in implicit form: t^2+ty+y^2 =7. This is a quadratic equation for y:
%    y^2+t\cdot y+(t^2-7)=0
% yielding y=\bigl(-t\pm\sqrt{t^2-4(t^2-7)}
% \bigr)/2
% =\bigl(-t\pm\sqrt{28-3t^2}
% \bigr)/2
% Plugging in t_0=1 gives 2=y_0=(-1\pm5)/2, hence we need "+":
% solution of IVP is y=\frac12\bigl(-t+\sqrt{28-3t^2}\bigr), interval of existence is \bigl(\,-\sqrt{\frac{28}3},
% \sqrt{\frac{28}3} \;\bigr)
%% Problem 1(d): t^2y+2ty^2+(t^2y+y^2)y'=0, y(0)=1
% This is an ODE M(t,y)dt+N(t,y)dy with M(t,y)=t^2y+2ty^2 and N(t,y)=t^2y+y^2
% Check whether this ODE is exact: need \partial_yM=\partial_t N:
%      \partial_yM = t^2+4ty,  \partial_t N=2ty
% Hence this ODE is not exact.
% Check whether integrating factor \rho(t) exists: we need that \frac{\partial_y M-\partial_tN}N is independent of y.
% Here we obtain \frac{\partial_y M-\partial_tN}N
% = \frac{t^2+2ty}{t^2y+y^2} which does depend on y. Hence no integrating factor \rho(t)exists.
% Check whether integrating factor \rho(y) exists: we need that -\frac{\partial_y M-\partial_tN}M is independent of t.
% Here we obtain -\frac{\partial_y M-\partial_tN}M
% = -\frac{t^2+2ty}{t^2y+2ty^2}=-\frac1y which does not depend on t. Hence we can find an integrating factor \rho(y) by solving \rho'(y)=-\frac{\partial_y M-\partial_tN}M \rho(y), i.e., \rho'(y)=-\frac1y\rho(y) yielding \rho^{-1}d\rho = -y^{-1}dy, \ln \rho = -\ln y, \rho(y)=y^{-1}.
% We now multiply the original ODE by the integrating factor \rho(y)=y^{-1} and obtain
%      (t^2 + 2ty)dt + (t^2+y)dy where we now have M(t,y)=t^2+2ty and N(t,y)=t^2+y
% We now check \partial_yM=\partial_t N:
%     \partial_yM = 2t ,    \partial_t N=2t  , hence this ODE is exact, and we can find H(t,y) with \partial_tH=M , \partial_yH=N.
% \partial_t H(t,y) = M =t^2+2ty gives H(t,y)=\frac{t^3}3+t^2y+g(y) where we still have to find the function g(y).
% \partial_yH(t,y)=N=t^2+y gives t^2+g'(y)=t^2+y \iff g'(y)=y (note that this does not depend on t), hence g(y)=\frac{y^2}2 (no need for constant, we just need to find one function H). We obtain
%      H(t,y)=H(t,y)=\frac{t^3}3+t^2y+\frac{y^2}2
% Note: we could also first use \partial_yH=N and then \partial_tH=M, this gives the same result. Just "guessing" H(t,y) and then showing that it works is not acceptable.
% The general solution in implicit form is H(t,y)=C. We plug in the initial condition t_0=0 and y_0=1 yielding C=\frac12.
% Hence the solution of the IVP in implicit form is \frac{t^3}3+t^2y+\frac{y^2}2 = \frac12.
% This is a quadratic equation for y:   y^2+2t^2\cdot y + \Bigl(\frac{2t^3}3-1\Bigr)=0, the quadratic formula gives
%    y(t) = \frac12 \left( -2t^2\pm\sqrt{4t^4-4\bigl(\frac{2t^3}3-1\bigr)}
% \right)
% =
%  -t^2\pm\sqrt{t^4-\frac{2t^3}3+1}
% Plugging in the initial condition t_0=0,\;y_0=1 gives 1=\pm\sqrt{1}, hence we need "+" for our initial condition.
% Hence the solution of the IVP is y(t)=-t^2+\sqrt{t^4-\frac{2t^3}3+1}. This solution exists as long as v(t):=t^4-\frac23 t^3+1\ge 0.
% For large |t| the term t^4 dominates, hence v(t)\to\infty as t\to\infty or t\to-\infty. We can find the local extrema of v(t) by solving v'(t)=0:   4t^3-2t^2=0 gives t^2(4t-2)=0, yielding t_1=0 and t_2=\frac12. At these points we have the function values v(0)=1 and v(\frac12)=\frac1{16}-\frac2{24}+1=1-\frac1{48}. Therefore we have v(t)\ge \frac{47}{48} for all t\in\mathbb{R}, hence the interval of existence is \mathbb{R}.
%% Problem 2(a)
% We have g(y)=y(y-1)(y-3)=0  for y=0, y=1, y=3. These are the stationary (constant) solutions of the ODE.

g = @(y) y*(y-1)*(y-3);
fplot(g,[-.4,3.2]); 
set(gca,'XAxisLoc','origin','YAxisLoc','origin'); box off    % draw axes thru origin
xlabel('y'); title('function g(y) = y(y-1)(y-3)')
f = @(t,y) g(y);
dirfield(f,0:.2:2,-.4:.2:3.4)
xlabel('t'); ylabel('y'); title('direction field for y''=y(y-1)(y-3)')

%% Problem 2(b)
% The function g(y) is zero for y_1^*=0,\; y_2^*=1,\; y_3^*=3. These are the stationary points.
% Note that the derivative g'(y) exists at the stationary points y_k^*. We can therefore use the theorem from class for \lim_{t\to\infty}y(t) and \lim_{t\to-\infty} y(t).
% For y_0=\frac12 we get \lim_{t\to-\infty} y(t) = 0,  \lim_{t\to\infty} y(t) = 1
% For y_0=2 we get \lim_{t\to-\infty} y(t) = 3,  \lim_{t\to\infty} y(t) = 1
% For y_0=4 we get \lim_{t\to-\infty} y(t)=3. For increasing t the solution y(t) will tend to \infty. There are two possibilities:
% Case 1: the solution exists for all t>0 and we have \lim_{t\to\infty}y(t)=\infty
% Case 2: the solution only exists for t<t_* and we have \lim_{t\to t_*}y(t)=\infty
% As explained in the hint on the ELMS page we have to check the improper integral \displaystyle I=\int_{y_0}^\infty \frac1{g(y)}dy. If this is infinite we have Case 1. If this is finite we have Case 2 with t_*=I.
% Here we get \displaystyle I = \int_4^\infty \frac1{y(y-1)(y-3)}dy \le \int_4^\infty \frac1{(y-3)^3}dy=
% \lim_{R\to\infty}\left[ \frac{-1}2(y-3)^{-2}\right]_4^R  is finite, hence we have Case 2.
% For y_0=-1 we get \lim_{t\to-\infty} y(t)=0. For increasing t the solution y(t) will tend to -\infty, this may either happen for t\to\infty (Case 1), or for t\to t_* (Case 2). For this we have to check the improper integral \displaystyle\int_{-\infty}^{-1}\frac1{g(y)}dy. With a similar argument as above we get that this integral is also finite, hence we have Case 2: the solution only exists for t<t_*, and \lim_{t\to t_*}y(t)=-\infty.
%
%% Problem 3(a)
% We are given the ODE (2xy^3+1)dx + (3x^2y^2+2y)dy=0 with M(x,y)=2xy^3+1 and N(x,y)=3x^2y^2+2y.
% Check whether the ODE is exact:
%       \partial_yM = 6xy^2 ,    \partial_xN=6xy^2
% We have \partial_yM = \partial_xN, hence the ODE is exact, and we can find H(x,y) with \partial_x H = M, \partial_yH=N.
% \partial_x H(x,y) = M = 2xy^3+1 gives H(x,y)=x^2y^3+x+g(y) where we still have to find the function g(y).
% \partial_yH(x,y)=N=3x^2y^2+2y gives 3x^2y^2+g'(y)=3x^2y^2+2y \iff g'(y)=2y (note that this does not depend on x), hence g(y)=y^2 (no need for constant, we just need to find one function H). We obtain
%      H(x,y)=x^2y^3+x+y^2
% The general solution in implicit form is H(x,y)=c with c\in\mathbb{R}.
%% Problem 3(b)

H = @(x,y) x^2*y^3+x+y^2;

% Plot the general solution in implicit form using fcontour:

fcontour(H,[-2 2 -2.5 1.5],'LevelStep',.4); colorbar

%% Problem 3(c): initial conditions (i) y(0)=1, (ii) y(0)=-1
% We plug in x_0=0,\;y_0=1

c1 = H(0,1)

% We plug in x_0=0,\;y_0=-1

c2 = H(0,-1)

% For both initial value problems we obtain the solution in implicit form H(x,y)=1. We can plot this using
% fcontour(H,[xl,xr yl yr],'LevelList',c):

fcontour(H,[-2 2 -2.5 1.5],'LevelList',1); hold on
plot(0,1,'bo'); 
plot(0,-1,'bo'); 
hold off; grid on; xlabel('x'); ylabel('y')

% We can see that the point (x,y)=(1,0) is on this curve: H(1,0)=1
% The gradient vector \bigl(\partial_x H,\partial_y H)=(M,N) at the point (1,0) is \bigl( M(1,0), N(1,0) \bigr) = \bigl(1,0\bigr) is horizontal. Since the gradient vector is orthogonal on the contour, the contour at this point is vertical.
% By looking at this graph:
% for initial condition y(0)=1 the solution seems to exist for t\in(-\infty,1)
% for initial condition y(0)=1 the solution seems to exist for t\in(-0.55,1)
%% Problem 3(d)
% We have M(x,y)+N(x,y)y'=0, hence
%          y'(x) = f(x,y(x))
% with f(x,y)=-M(x,y)/N(x,y).

f = @(x,y) -(2*x*y^3+1)/(3*x^2*y^2+2*y);
dirfield(f,-2:.2:2,-2.5:.2:1.5); 
hold on; xlabel('x'); ylabel('y'); title('(2xy^3+1)+(3x^2y^2+2y)y''=0, initial conditions y(0)=1 and y(0)=-1')

%% First initial condition: y(0)=1

x0 = 0
y0 = 1;
plot(x0,y0,'bo')
[xs,ys] = ode15s(f,[0,2],y0);    % solve from x=0 to 2
[xs(end),ys(end)]

% The solution ends close to x\approx 1 and y\approx 0. Note that ode15s approximates the solution, so the computed values are not exact.

plot(xs,ys,'b')
[xs,ys] = ode15s(f,[0,-2],y0);   % solve from x=0 to -2
plot(xs,ys,'b')

%% Second initial condition: y(0)=-1

y0 = -1;
plot(x0,y0,'go')
[xs,ys] = ode15s(f,[0,2],y0);    % solve from x=0 to 2
[xs(end),ys(end)]

% The solution ends close to x\approx 1 and y\approx 0.

plot(xs,ys,'g')
[xs,ys] = ode15s(f,[0,-2],y0);   % solve from x=0 to -2
[xs(end),ys(end)]

% The solution ends close to x\approx -0.55 and y\approx -2.18.

plot(xs,ys,'g')
hold off