How can I find an approximation for y(1) using MATLAB without ODE solvers? (Euler's Method & Matrices) Task:
Given \begin{align*} y(t+h)&\approx \underbrace{\left(\begin{matrix}1-h & 5h & h\\3h & 1-h & 0\\0 & -th & 1+h\end{matrix}\right)}_{F(t,h)}y(t)+\underbrace{\left(\begin{matrix}0\\ht^2\\0\end{matrix}\right)}_{g(t,h)} \end{align*}
and $h=0.1$
and $y(0)={\left(\begin{matrix}0\\0\\0\end{matrix}\right)}$
Question:
How can I find an approximation for $y(1)$ using MATLAB, without using ODE solvers (i.e. write an algorithm)? 
Extra Info. that might be useful:
$y(t)={\left(\begin{matrix}y_1(t)\\y_2(t)\\y_3(t)\end{matrix}\right)}$
$y'(t)={\left(\begin{matrix}-1 & 5 & 1\\3 & -1 & 0\\0 & -t & 1\end{matrix}\right)}y(t)+{\left(\begin{matrix}0\\t^2\\0\end{matrix}\right)}$

You have written down one step of the Euler method. Now apply this 10 times to get from 0 to 1.
At $t=0$ you compute $y(0.1)$ with $h=0.1$, from that you compute $y(0.2)$ etc.
Be advised that the Euler method has a global error (more relative than absolute) of $O(h)$, thus with $h=0.1$ you can expect an error of magnitude of $10\%$, $y_h(1)-y(1)=C·h+O(h^2)$. To get a better measure of the error, perform the same computation with $h=0.2$ and $h=0.05$. You should get that 
\begin{align}
y_{0.2}(1)&-y_{0.1}(1)&&\approx &2·(y_{0.1}(1)&-y_{0.05}(1))\\
&\approx &&&&\approx \\
C·0.1&+O(0.01)&&\approx& 2·(C·0.05&+O(0.0025))
\end{align}
and if that is true then a better guess for the exact value is obtained from
$$
2·y_h(1)-y_{2h}(1)=y(1)+O(h^2).
$$
See Richardson extrapolation.
