Uwe,
In my testing (using octave since I don't have matlab, but I hope it's similar), using = :exports none :results output raw= seems to produce the desired output:
#+begin_src octave :exports results :results output raw
close all
N = 3; % number of chebyshev nodes
n = 1; % polytropic index
iters = 1; % iterations
j=[0 1 2 3];
z=cos((pi*j)/3);
y=1-(0.5*(z+1)).^2;
a=n*(8);
a2=eye(3+1)*a;
disp('\begin{align*}')
disp('A_{1,k-1}&=')
disp('\begin{pmatrix}')
fprintf('%g &%g &%g &%g\\\\\n',a2')
disp('\end{pmatrix}=OK\\')
disp('\end{align*}')
#+end_src
#+RESULTS:
\begin{align*}
A_{1,k-1}&=
\begin{pmatrix}
8 &0 &0 &0\\
0 &8 &0 &0\\
0 &0 &8 &0\\
0 &0 &0 &8\\
\end{pmatrix}=OK\\
\end{align*}
Which renders the matrix correctly. I first tried with =:exports results=, but that for some reason results in the results being included twice in the LaTeX output. I'm not sure why.
Hope this helps,
--Diego