Consider the following equation: $u_t=div\Bigg(\begin{pmatrix}a& b\\b & c\end{pmatrix}\nabla u\Bigg)$. Simplifying we get $u_t=au_{xx}+2bu_{xy}+cu_{yy}$. I am trying to apply this equation on Images and see the changes. It is clear that if the matrix is identity matrix then we get back the heat equation. I have written a matlab code for this. There is some problem with the code. For values of $a,b,c < 1$, it works fine. But for values $>1$ I get a weird output. I am attaching the matlab code and the output here:
% To apply affine diffusion equation u_t=div(A*grad u) on images
% where A is a positive definite matrix
% Accept Image from the user
disp('Enter Image')
im=uigetfile('*.jpg','Select the Image file');
% input image file and store in u1,u2,u3
u1=imread(im);
u1=double(u1);
% Accept the matrix from the user
a=input('Enter a: ');
b=input('Enter b: ');
c=input('Enter c: ');
T1=input('Enter Stop Time: ');
% Display Original Image
subplot(1,2,1), imshow(uint8(u1))
title('Original image');
% store size in m and n
[m,n,~]=size(u1);
% step size along t
dt=0.25;
for t = 0:dt:T1
% finite difference approximation for u_xx
u_xx = u1(:,[2:n n],:) - 2*u1 + u1(:,[1 1:n-1],:);
% finite difference approximation for u_xy
u_xy = ( u1([2:m,m],[2:n,n],:) + u1([1,1:m-1],[1,1:n-1],:) - ...
u1([1,1:m-1],[2:n,n],:) - u1([2:m,m],[1,1:n-1],:) ) / 4;
% finite difference approximation for u_yy
u_yy = u1([2:m m],:,:) - 2*u1 + u1([1 1:m-1],:,:);
% Calculate RHS
div_mat_u=a*u_xx+2*b*u_xy+c*u_yy;
u1= u1 + dt*(div_mat_u);
temp=u1;
end
subplot(1,2,2), imshow(uint8(temp))
str1=sprintf('T=%d',T1);
title(str1);
clc;
And the sample output with the matrix \begin{pmatrix}2& 0\\0 & 2\end{pmatrix} with stopping time $T=10$ is here:
The sample output with the matrix \begin{pmatrix}0.9& 0.2\\0.2 & 0.8\end{pmatrix} with stopping time $T=10$ is 
What am I doing wrong here ?
You can do a von Neumann analysis to check stability of the scheme. I will illustrate on the 1D heat equation
$$u_t = au_{xx}.$$
The explicit scheme is
$$u^{n+1}_j = (1-2s)u^n_j + s(u^n_{j+1}+u^n_{j-1}),$$
where $s=adt/dx^2.$ The basics of a von Neuman analysis are to take $u^n_j = e^{ijk}$ (here, $i=\sqrt{-1}$) and look for the amplification factor $\lambda$ for which $u^{n+1}_j = \lambda u^n_j$. The scheme is stable if and only if $|\lambda |\leq 1$ for all $k$. At a deeper level, we are expanding the solution in a Fourier basis.
Plugging this in and cancelling the $e^{ijk}$ from both sides we get
$$\lambda = (1-2s) + s(e^{ik}+e^{-ik}) = (1-2s) + 2s\cos(k).$$
We just have to look at the extremes, $\cos(k)=1$ and $\cos(k)=-1$ to get $$1-4s \leq \lambda \leq 1.$$ Thus we need $1-4s \geq -1$ to ensure $-1 \leq \lambda \leq 1$. This is equivalent to $s\leq 1/2$ or $dt \leq dx^2/(2a)$. Hence, if you increase $a$ you need to decrease $dt$ for stability.
A similar analysis can be carried out in 2D, but it is a bit more tedious. In your code, you may want to normalize the matrix $A$ given by the user (say, divide $A$ by a constant so that the largest entry is 1. Actually, probably the trace of $A$ should play the role of $a$ from the analysis above). Then you can keep the same time step. You also need to ensure $A$ is positive definite or it is not a heat equation (it becomes a backwards heat equation and will be ill-posed for other reasons).
As another remark, if you are interested in image processing, you should look into nonlinear diffusion equations (such as the Perona-Malik anisotropic diffusion). The idea is to make $A$ dependent on the image so that smoothing occurs along edges and not across edges (thereby preserving edges).
EDIT: Let me add a short discussion about implicit schemes. Implicit schemes discretize $u_{xx}$ at the next time step $n+1$, hence the scheme is
$$u^{n+1}_j = u^n_j + s(u^{n+1}_{j+1}-2u^{n+1}_j+u^{n+1}_{j-1}).$$
Rearranging we have
$$(1+2s)u^{n+1}_j = u^n_j + s(u^{n+1}_{j+1}+u^{n+1}_{j-1}).$$
Running the same von Neumann analysis you get
$$(1+2s)\lambda = 1 + 2s\lambda \cos(k).$$
Hence
$$\lambda = \frac{1}{1+2s - 2s\cos(k)}.$$
The denominator is always $\geq 1$ hence $0\leq \lambda \leq 1$ independent of $s$. Hence the scheme is unconditionally stable.
You can also use a similar implicit scheme for your PDE. This allows you to take arbitrarily large time steps (the only concern is accuracy), but each iteration is more expensive comptuationally since it involves solving a linear system (it's not explicit).