Here is a sample code, that recovers $\Phi_2$ and $\Phi_4$ from

https://math.mit.edu/~drew/ClassicalModPolys.html

using sage code. (Note that the polynomial displayed in the above link is $\Phi_2$, not $\Phi_4$ as stated.)
One needs the **kohel** database. (In manjaro i installed it from the AUR, the name of the package was `sage-data-kohel`

, referenced here, sage-data-kohel so i can not state here certainly that `sage -i database_kohel`

does the job. A similar command, possibly with a suffix should install the database properly.)

Then

```
[dan@k7 ~]$ sage
┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 8.1, Release Date: 2017-12-07 │
│ Type "notebook()" for the browser-based notebook interface. │
│ Type "help()" for help. │
└────────────────────────────────────────────────────────────────────┘
sage: PHI = ClassicalModularPolynomialDatabase()
sage: var( 'X,Y' );
sage: PHI[2]
-j0^2*j1^2 + j0^3 + 1488*j0^2*j1 + 1488*j0*j1^2 + j1^3 - 162000*j0^2 + 40773375*j0*j1 - 162000*j1^2 + 8748000000*j0 + 8748000000*j1 - 157464000000000
sage: PHI[2](X,Y)
-X^2*Y^2 + X^3 + 1488*X^2*Y + 1488*X*Y^2 + Y^3 - 162000*X^2 + 40773375*X*Y - 162000*Y^2 + 8748000000*X + 8748000000*Y - 157464000000000
```

reproduces the (symmetrical) $\Phi_2(X,Y)$ polynomial from the above link. To have the data from the links

we can then ask for

```
PHI = ClassicalModularPolynomialDatabase()
for N in (2, 3, 4):
print "N = %s" % N
data = [ ( mono.degrees(), coeff ) for coeff, mono in PHI[N] ]
data . sort()
for degrees, coeff in data:
# we show only the half of the coefficients of the symmetrical pol PHI[N]
if degrees[0] >= degrees[1]:
print list(degrees), coeff
print
```

and the results are:

```
N = 2
[0, 0] -157464000000000
[1, 0] 8748000000
[1, 1] 40773375
[2, 0] -162000
[2, 1] 1488
[2, 2] -1
[3, 0] 1
N = 3
[1, 0] 1855425871872000000000
[1, 1] -770845966336000000
[2, 0] 452984832000000
[2, 1] 8900222976000
[2, 2] 2587918086
[3, 0] 36864000
[3, 1] -1069956
[3, 2] 2232
[3, 3] -1
[4, 0] 1
N = 4
[0, 0] 280949374722195372109640625000000000000
[1, 0] -364936327796757658404375000000000000
[1, 1] -94266583063223403127324218750000
[2, 0] 158010236947953767724187500000000
[2, 1] 188656639464998455284287109375
[2, 2] 26402314839969410496000000
[3, 0] -22805180351548032195000000000
[3, 1] 12519806366846423598750000
[3, 2] -914362550706103200000
[3, 3] 2729942049541120
[4, 0] 24125474716854750000
[4, 1] 1194227244109980000
[4, 2] 1425220456750080
[4, 3] 80967606480
[4, 4] 7440
[5, 0] -8507430000
[5, 1] 561444609
[5, 2] -2533680
[5, 3] 2976
[5, 4] -1
[6, 0] 1
```

The relevant py-file may be located as follows:

```
sage: M = sage.databases.db_modular_polynomials
sage: M.__file__
'/usr/lib/python2.7/site-packages/sage/databases/db_modular_polynomials.pyc'
```

and the py (not its compiled version, with suffix pyc) starts with:

```
"""
Database of Modular Polynomials
"""
#######################################################################
# Copyright (C) 2006 William Stein <wstein@gmail.com>
# Copyright (C) 2006 David Kohel <kohel@maths.usyd.edu.au>
# Copyright (C) 2016 Vincent Delecroix <vincent.delecroix@labri.fr>
#
# Distributed under the terms of the GNU General Public License (GPL)
# as published by the Free Software Foundation; either version 2 of
# the License, or (at your option) any later version.
# http://www.gnu.org/licenses/
#*****************************************************************************
from __future__ import print_function, absolute_import
def _dbz_to_string(name):
```

and so on. Hope, this is a good starting point to get the modular polynomials, the canonical equations for the modular curves in the family $X_0(N)$. (For me, `PHI[19]`

was printed, `PHI[20]`

was missing, `LookupError: filename /usr/share/kohel/PolMod/Cls/pol.020.dbz does not exist`

.)