I dont have a background in maths and this is my first time dealing with pdes, trying to solve $$ \begin{matrix} \frac{\partial n_c}{dt} \\ \frac{\partial n_p}{dt} \\ \end{matrix} =\begin{pmatrix} D_c \nabla^2-k_c-k & 0 \\ k & D_c \nabla^2-k_p \\ \end{pmatrix} \begin{pmatrix} n_c(r,t) \\ n_p(r,t) \\ \end{pmatrix} $$
After some attempt I could see that $n_c$ has an analytical solution, but I dont think $n_p$ has one. What are the tools that I could use to tackle this problem? I'm proficient with Python but could also try learning other language for solving.
Assuming $k, k_c, k_p, D_c$ are constants, you can write down an analytical solution for $n_p$ reasonably easily. As you have observed, the PDE for $n_c$ is decoupled from the PDE for $n_p$, and boils down to the equation $$ \frac{\partial n_c}{\partial t} = (D_c\nabla^2 - k_c - k)n_c $$ which is effectively a heat equation and can be solved analytically with a variety of techniques (such as separation of variables/Fourier analysis) given appropriate initial/boundary conditions.
Once you have solved for $n_c$, you can now solve for $n_p$, because the equation for $n_p$ is of the form $$ \frac{\partial n_p}{\partial t} = kf(r,t) + (D_c\nabla^2 - k_p)n_p $$ where $f(r,t) = n_c(r,t)$. This is an inhomogeneous linear PDE, and assuming you have reasonable initial conditions, it can be solved using Duhamel's principle. Note that you can use the explicit analytical solution for $n_c$ in the PDE for $n_p$, so there is no issue arising from the equations being coupled: you solve for $n_c$ first, then plug it back into the PDE for $n_p$ and solve from there.