Question: I have created a large positive definite (thus symmetric) matrix in Matlab. For sake of simplicity, let's assume that the matrix has the following form:

$K= \left[ \begin{array}{ccc} k\_{11} & \cdots & k\_{1n} \\ \vdots & \ddots & \vdots \\
k\_{n1} & \cdots & k\_{nn}\end{array} \right]$.

Now I want to raise it to a power in an elementwise way, let's say to the power of $\lambda$.

Since $K$ is known, while $\lambda$ is unknown and has to be estimated via an optimization routine, I save beforehand the matrix $K$. In the optimization routine I have to compute at each iteration the value of the matrix raised to the power of $\lambda$.

My question is the following: Is there a way to perform this operation in a fast way, by exploiting the properties of the matrix $K$?

What I have done until now:

1st approach: Compute $K^\lambda$ in a straightforward way, that is to say, using command "K.^\lambda" in matlab. 

2nd approach: Use two for loops, where the first one iterates from i=1:n and the second one from j=i:n. This will create an upper diagonal matrix and due to symmetry, I can use triu command to get back the full matrix.

Nonetheless, both approaches are rather time consuming, due to the large size of the matrix. Does anyone have any better idea about how I could perform this operation?

Answer: 
It is entirely possible, although I cannot tell by looking at it, that your system becomes stiff near $\eta=0$. I say this because Mathematica's NDSolve integrates this up to $\eta=3.4\times 10^{-9}$, then complains about a singularity or stiffness. 

The midpoint method is not A-stable, so that may be the source of inaccuracy: it is unsuitable for stiff problems. The ODE has some kind of singularity near $0$, so maybe the solution blows up because it hits the singularity due to inaccuracy before the method's instability becomes apparent. What do you get if you use a larger time step, like $0.01$ or $0.001$?

Try an A-stable method, like the trapezoid rule:
$$ y\_{n+1} = y\_n + \frac{h}{2}(F(x\_n,y\_n) + F(x\_{n+1}, y\_{n+1})). $$
You'd need to solve the implicit system of nonlinear equations here.

Another approach, since you know the asymptotic behaviour of the ODE solution, is to write the unknown functions as perturbations around the known asymptotic behaviour (e.g., $\frac1\gamma\eta+\eta^2\hat\phi(\eta))$, etc., as in a Taylor series), which might make the resulting ODE in $\hat\phi$ easier to solve by (maybe) avoiding the numerical errors in the denominator.
My suggestion is to analytically detect the singularities and then to switch form explicit solver to implicit one around the singularity.
From this equation
$
\psi ' = \dfrac{\phi ' + 2 \frac{\phi}{\eta}}{\eta-\phi}\psi.
$
it follows that the singularity is around $\eta-\phi$. Therefore, you should require that 
$
\phi ' = - 2 \frac{\phi}{\eta} 
$
with a solution $\phi=\frac{c}{{\eta}^{2}}$.
Therefore 
$
\eta = \frac{c}{{\eta}^{2}}
$
and the critical values is $\eta^\*= ^3\sqrt{c}$.

The second equation
$
\phi'=\dfrac{\frac{1}{\gamma} f'}{\psi(\eta-\phi)}-\dfrac{3}{2}\dfrac{\phi}{\eta-\phi}
$
will be then
$
 f' = \dfrac{3 \gamma \psi}{2} \phi 
$

The first equation
$
f'=\dfrac{f[-3\eta + (3+\gamma/2)\phi - 2\gamma \dfrac{\phi ^2}{\eta}]}{(\eta-\phi)^2-f/\psi}
$
will become
$
f'=\psi(-3\eta + (3+\gamma/2)\phi - 2\gamma \dfrac{\phi ^2}{\eta}) 
$