Solving Differential Equations theoretically and using matlab i am trying to solve the initial value and elliptic boundary value problems below. but now i need some help solving them using matlab. for the elliptic problem, any method is ok, but for the initial value problem, i need to use finite difference method (any choice). any help with code, or guidance to functions is greatly appreciated in advance. i have go Matlab R2011b. below are the equations:
1.
\begin{align}
&y''(t) = exp(-t) + 3y(t), \quad 0 \leq t \leq 2, \\
&y(0) = 1, \\
&y'(0) = -1
\end{align}
by finite diff method of order $O(h^4)$
2.
\begin{align}
&-U_{yy}(y,x) - U_{xx}(y,x) = 2(\pi^2) \sin(\pi x) \cos(\pi y), \quad y,x \in [0,1], \\
& U_y(0,x) = U_y(1,x) = 0, \\
&U(y,0) = U(y,1) = 0. 
\end{align}
Any help is greatly appreciated. a pointer to function, or similar problem is as well appreciated. thanks.
 A: I will describe in summary form how I would go on solving these problems:
Problem 1:
First note that you can re-formulate your initial value problem (IVP) as follows:
$$ 
Y'(t) = F(t,Y(t)), \quad  t > 0,
$$ where $Y(t) = (w(t),y(t)), $ and $w(t) := y'(t)$ has been defined for convenience. Then, the initial conditions are given as: $Y(0) = (w(0),y(0)) = (y'_0 , y_0)$, which are known values. If a 4th order accurate method is required, then you will need at least three known points of the solution, at least, i.e., the numerical solution for the points $(t_0,t_1,t_2)$. You can proceed according to many ways. One of them would be: 


*

*Solve the first three points using an Euler's scheme:
$$Y_{n+1} - Y_n = F(t_{n+s},Y_{n+s}), \quad n = {1,2}, \quad Y_0 = (y'_0,y_0),$$ where $s=0$ if the scheme is explicit and $s=1$ if implicit. 

*Then use a BDF-4 method to march in time:
$$ Y_{n+4} - \frac{12}{25} h \, F(t_{n+4} , Y_{n+4}) = \frac{48}{25} Y_{n+3} - \frac{36}{25} Y_{n+2} + \frac{16}{25} Y_{n+1} - \frac{3}{25} Y_n, \quad n = {0, 1, \ldots},$$ and you are done. You can see here how to solve a given 2nd order IVP in Matlab using a forward Euler's scheme. More information about Matlab's solvers can be found here, but I would recommend to implement your own codes. Other suitable methods for solving this problems can be: Runge-Kutta, Milne's method, Adams-Bashforth, Adams-Moulton or other predictor-corrector formulae.
Problem 2:
Your problem may be rewritten as follows: 
$$ - \nabla^2 u = F(x,y), \quad (x,y) \in [0,1]\times[0,1], $$ with the correspondent Dirichlet conditions in $x$ and Neumann's in $y$. A finite difference method can be developed and implemented in Matlab considering the uniform mesh:
$$ \mathcal{M}  = \{ (x_i,y_j) \in \mathbb{R}^2 \ x_i = x_0 + i \Delta_x , \ y_j  = y_0 + j \Delta_y  \}, $$ where $x_0 = y_0 = 0,$ $\Delta_x = X/N_x, \ \Delta_y = Y/N_y$ with $X$ and $Y$ the lengths of the $x_i$-intervals and $N_{x_i} + 1$ the number of nodes in which the $x_i$ coordinates are discretized. Then, if we assume $U_{ij} \approx u(x_i,y_j)$, the equation can be formulated as follows:
$$ - \frac{1}{\Delta x^2} (U_{i+1,j} - 2 U_{ij} + U_{i-1,j} ) - \frac{1}{\Delta y^2}(  U_{i,j+1} - 2 U_{ij} + U_{i,j-1}  ) = F_{ij},  $$ with $F_{ij} = F(x_i,y_j)$. A second order centered finite difference scheme has been used to discretize the Laplace's operator. The boundary conditions are to be also discretized: 
$$\begin{align}
U_{0,j} & = 0, \\
U_{N_x,j} & = 0, \\
U_{i,-1}  & = U_{i,1}, \\
U_{i,N_y+1} & = U_{i,N_y-1}, \\
\end{align} $$ note that the two last points are notes out of the boundary and must be substitudes in the original equation. Once you write the FDM equation for the corresponding indices of the boundary conditions, you will end up with a system of equations which can be written of the form:
$$ \alpha U_{i-1,j} + \beta U_{ij} + \gamma U_{i+1,j} + \delta_1 U_{i,j-1} + \delta_2 U_{i,j+1}= \eta F_{ij}, $$ where some $\alpha, \beta,$ etc. coefficients may vary when applying the Neumann boundary conditions. This system admits a matrix form as follows:
$$ \mathbf{A} \mathbf{U} = \mathbf{F}, \quad \mathbf{A} = \left(
\begin{array}{ccccc}
 A & D_2 & \cdots  & O & O \\
 D_1 & A & \cdots  & O & O \\
 \vdots  & \vdots  & \ddots & \vdots  & \vdots  \\
 O & O & \cdots  & A & D_2 \\
 O & O & \cdots  & D_1 & A
\end{array}
\right), 
$$ 
$$\mathbf{U} = (U_{10},U_{20}, \ldots, U_{N_x 0} | \ldots |U_{11},U_{21}, \ldots, U_{N_x1} | \ldots | U_{1N_y},U_{2N_y}, \ldots, U_{N_x-1, N_y} )^T,
$$ $A$ is a tridiagonal matrix with $\alpha,\beta,\gamma$ coefficients and $D_k$ are diagonal matrices with $\delta_k$ coefficients. This system is sparse, so it can be solved in a cleaver way using some sparse solver with the help of Matlab. 
Edit: 
Here and here you can find two examples (of the same author) for solving the Laplace's equation with Matlab. The question is: how to assemble the big matrix $\mathbf{A}$? Remember that this matrix is formed in itself by the smaller matrices $A$ and $D_{1,2}$, which are $(N_x - 1)$-squared matrices so the total system becomes a $(N_x-1) \times (N_y +1)$-squared set. I'll let you to tell me why these are the sizes of the correspondent matrices (hint: look at the BCs). 
I like to think that $i$ are the local indices while $j$ may be defined as global indices. Being that the case, you readily follow that every magnitudes with subscripts $j+1$ and $j-1$ are superdiagonal matrices $D_2$ and subdiagonal matrices $D_1$ which are global in the sense that they are full matrices The coefficients of $D_k$ are diagonal since they are subscripted with $i$. Now, as the other coefficients have the subscript $j$, it means that they remain global diagonal matrices, $A$, which coefficients are given by the local indices $i-1, \ i $ and $i+1$: $\alpha, \ \beta$ and $\gamma$.
Hope this helps!
If you give me some minutes, I could add some hints for solving both problems analytically. 
Cheers!
