def maxTime():
    def heatmap(ix,iy,iz,borders):
        levels = MaxNLocator(nbins=50).tick_values(iz.min(), iz.max())
        # pick the desired colormap, sensible levels, and define a normalization
        # instance which takes data values and translates those into levels.
        cmap = plt.get_cmap('Greys')
        norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

        fig, (ax0, ax1) = plt.subplots(nrows=2)

        im = ax0.pcolormesh(ix,iy, iz, cmap=cmap, norm=norm)
        fig.colorbar(im, ax=ax0)
        ax0.set_title('pcolormesh with levels')


        # contours are *point* based plots, so convert our bound into point
        # centers
        cf = ax1.contourf(ix[:-1, :-1] + 0.01/2.,
                          iy[:-1, :-1] + 0.01/2., iz, levels = borders,
                          cmap=cmap)
        fig.colorbar(cf, ax=ax1)
        ax1.set_title('contourf with levels')

        # adjust spacing between subplots so `ax1` title and `ax0` tick labels
        # don't overlap
        fig.tight_layout()
        plt.show()


    x = np.arange(0.1, 1, 0.1)
    y = np.arange(0.1, 1, 0.1)

    t = sp.Symbol('t', positive=True)
    p = sp.symbols('p')
    Ki = sp.symbols('Ki', real = True)
    Kp = sp.symbols('Kp', real=True)
    a = sp.symbols('a', real=True)
    b = sp.symbols('b', real=True)
    c = sp.symbols('c', real=True)
    expression = (p*Kp+Ki)*a/(p*(b*c*p+a*a)+(p*Kp+Ki)*a)

    orig1 = sp.inverse_laplace_transform(expression, p, t)
    orig2 = sp.inverse_laplace_transform(expression/p, p, t)
    x, y = np.meshgrid(x, y)
    tp = np.empty([9,9])
    mp = np.empty([9,9])
    for i in range(9):
        for j in range(9):
            tp[i][j] = nsolve(simplify(orig1.subs([(Ki,x[i][j]),(Kp,y[i][j]),(a,0.03),(b,2.44),(c,0.00003)])),t,(0,1))
            print(orig2.subs([(Ki,x[i][j]),(Kp,y[i][j]),(a,0.03),(b,2.44),(c,0.00003),(t,tp[i][j])]))

    tp = tp[:-1, :-1]
    heatmap(x,y,tp,[0,15,500])




maxTime()