In order to simulate the following equation using FDM $$u_t(t,x)-u_{xx}(t,x)=0, \quad (t,x) \in (0,1)\times (0,1)$$ $$(u_t(t,x)-u_{x}(t,x))\rvert_{x=0}=0, \quad t \in (0,1)$$ $$(u_t(t,x)+u_{x}(t,x))\rvert_{x=1}=0, \quad t \in (0,1)$$ $$u(0,x)=u_0, \quad x \in (0,1)$$ $$u(0,0) =\mu_0, \quad u(0,1)=\mu_1,$$ where $f,g,u_0$ are given functions and $u_1,u_2$ are constant. I'm struggling in the discretization to obtain the corresponding linear system.
For the first equation taking $u_j^n=u(t_n,x_j)$ and using explicit Euler scheme, I found $$u_j^{n+1}=ru_{j-1}^n+(1-2r)u_j^n +ru_{j+1}^n,$$ with $r=\dfrac{\Delta t}{(\Delta x)^2}$, but I don't know how to discretize the second and third equations. For the $x$-derivatives I think I should use centered formula by adding two ghost points as in Neumann b.c, and the same for time derivatives but this is a bit complicated and yields much of equations.
Maybe I should treat the equations like a coupled system through the boundary but I don't know how to do that. What is the right FDM discretization. Thanks for any help.
PS: My background to FDM is very simple (just theoretical) and this is my first time to simulate PDE. So, what is your advice to this in a right way? (methods, languages, programs...etc).