Using python, can you step by step break down how to code/model the following pde's and then plot

a/t = D_a*a + _aF(a,h) _a*a

h/t = D_h*h + _hF(a,h) _h*h

D_a = 0.005, D_h = 0.2, _a = 0.005 _h = 0.02, _a = 0.005,

_h =0.015, K_h = 0.01, and = 0.1

F(a,h) = a^2/(Kh + h)(1+ a^2)

a=Daa+aF(a,h)aa,h=Dhh+hF(a,h)hh. F(a,h)=(Kh+h)(1+a2)a2 As parameter values we chose Da=0.005,Dh=0.2,a=0.005h=0.02,a=0.005,h= 0.015,Kh=0.01, and =0.1. Here it should be noted that the actual parameter values are arhitrary and not based on experimentally measured quantities. However, the system is robust