How does one compute the space of algebraic global differential forms $\Omega^i(X)$ on an affine complex scheme $X$? In 1963 Grothendieck introduced the algebraic de Rham cohomolog in a letter to Atiyah, later  published in the Publications Mathématiques de l'IHES, N°29.
If $X$ is an algebraic scheme over $\mathbb C$ the definition is via the hypercohomology of the complex of Kähler differntial forms: $$H^i_{dR}(X)=\mathbb H^i(X,\Omega_{X/\mathbb C}^{\bullet})$$ If in particular$X$ is affine: $$H^i_{dR}(X)=h^i(\Omega_{X/\mathbb C}^\bullet (X))$$ If moreover $X$ is affine and smooth: $$H^i_{dR}(X)=H^i_{sing}(X(\mathbb C),\mathbb C)$$ so that we get a purely algebraic means of computing the cohomology of the set of closed  points $X(\mathbb C)$ of our scheme  provided with its transcendental topology (derived from the topology of the metric space $\mathbb C$).
In order to understand this algebraic de Rham cohomology  I tried to compute some examples in the non smooth case, in particular  for the cusp $V\subset \mathbb A^2_\mathbb C$ given by  $y^2=x^3$.
The slightly surprising result is that $H^1_{dR}(V)$ is the one-dimensional complex vector space  with basis the class of the differential form $x^2ydx$.
This is in stark contrast with the fact that $H^1_{sing}(V(\mathbb C),\mathbb C)=0$, since $V(\mathbb C)$ is homeomorphic to $\mathbb C$. 
I made some very stupid preliminary mistakes: in particular I hadn't realised   at the start of my calculation that $d(x^2ydx)=x^2dy\wedge dx=0$, as follows from the equality $2ydy=3x^2dx$ (and  wedge-multiplying by $dy$).
In other words $\Omega^2(V)$ was more complicated than I thought: although I knew  the formula for, say, the differentials of order one for a hypersurface $X=V(f)\subset \mathbb A^n_\mathbb C $, namely $\Omega^1(X)=\frac { \sum \mathcal O(X)dx_j}{\langle \sum \frac {\partial f}{\partial x_j} dx_j\rangle }$, I realized that I didn't know the algorithm for computing higher order global differentials.
This is the motivation for my
ACTUAL QUESTION
Given the  subscheme $$X=V(f_1(x),f_2(x),\cdots, f_r(x))\subset \mathbb A^n_\mathbb C \quad (f_j(x)\in \mathbb C[x_1,\cdots, x_n]),$$  how does one compute the space of algebraic global differential forms $\Omega^i_{X/\mathbb C}(X)$ ?
(Ideally I'd love to learn an algorithm which for example would calculate $\Omega^2(V)$ for the cusp above in a picosecond)
Edit
The software showed me a related question that my computation for the cusp solves.
I have used this opportunity to write down the details of that computation as an answer to the related question.
 A: This is not an answer, but is too long for a comment.
Macaulay 2 can compute the exterior powers of the cotangent sheaf of a projective variety, and it can compute the de Rham cohomology of the complement of an affine hypersurface.  I link to the M2 book for a more complete description of the first approach, and to the M2 documentation for the de Rham cohomology of a complement.  Although I think both approaches can be used to answer your question completely, I have not been able to figure it how to do it either way.
M2 uses the conormal sequence to find a presentation for the cotangent sheaf.  Concretely, this is the Jacobian matrix of the homogeneous forms that cut out your ideal.  For example, to compute the sheaf $\Omega^2(X)$ where $X$ is the subscheme defined by the homogeneous polynomial $y^2z-x^3$, one uses the following command:
X=Proj(QQ[x,y,z]/(y^2*z-x^3));
omega2=cotangentSheaf(2,X)

which returns the presentation matrix
cokernel {0} | 0 y 0 0 x2 0  -3x2 0   0   -3yz z4  0   |
         {1} | 0 0 y 0 0  x2 -2y  0   0   -2x  0   z4  |
         {3} | 0 0 0 0 0  0  0    -2y 3x2 0    -2x 3yz |

                                   1       1           1
coherent sheaf on X, quotient of OO   ++ OO  (-1) ++ OO  (-3)
                                   X       X           X

I do not know how to build the exterior derivative to relate the various wedge powers of omega.  However, once that's done, M2 can do the hypercohomology spectral sequence.  Maybe someone who knows how to do this will answer.
The second approach can compute something similar to the de Rham cohomology you want, since the cohomology of a complement is related by a long exact sequence to the Borel-Moore homology, which is the relative homology of the one-point compactification of $V$ rel the added point.
Here is how to compute the cohomology of the complement of $y^2=x^3$ in $\mathbb{C}^2$:
loadPackage("DModules");
R = QQ[x,y];
f = y^2-x^3;
deRham f

which returns
                 1
HashTable{0 => QQ }
                 1
          1 => QQ
          2 => 0

You can get more detailed information about the differential forms representing these classes by using
deRhamAll f

I hope someone else can explain the proper D-module calculation to answer the question!
A: $\require{AMScd}$
$\require{AMStex}$
$\newcommand{\CC}{{\mathbb C}}$
$\newcommand{\PSP}{\mathbf{\mathbb P}}$
$\newcommand{\AFF}{\mathbf{\mathbb A}}$
$\newcommand{\Ohol}{{\mathcal O}}$
$\newcommand{\Hdr}{H_{\text{DR}}}$
$\newcommand{\HH}{{\mathbb H}}$
$\newcommand{\cover}[1]{{\mathfrak #1}}$
$\newcommand{\Cc}{\check{C}}$
$\newcommand{\codim}{\text{codim}}$
The following remarks outline a strategy for computing $\Hdr^q(X)$ for every quasiprojective variety $X \subseteq \PSP^n_\CC$. I call it a strategy rather than an algorithm, as I think some parts of the calculation still need a lot of thought, to be implemented as an algorithm, especially where spectral sequences appear:
Let $X/k$ be an algebraic variety over $k = \CC$. Consider the sequence
$(\Omega^\bullet_{X|k})$
$$
\Ohol_X = \Omega^0_{X|k} \xrightarrow{d^0} \Omega^1_{X|k} \xrightarrow{d^1} \Omega^2_{X|k} 
\xrightarrow{d^2} \cdots
$$
where $\Omega^p_{X|k} = \bigwedge^p \Omega_{X|k}$ and $d^p:\Omega^p_{X|k} \to \Omega^{p+1}_{X|k}$
stands for the exterior derivative. It is $k$--linear but not $\Ohol_X$--linear.
One defines the $i$--th de Rham cohomology of $X$ as
$$
\Hdr^i(X) = \HH^i(\Omega^\bullet_{X|k})
$$
where $\HH^i(L^\bullet)$ is the $i$--th hypercohomology of the complex $L^\bullet$.
If $X$ is an affine variety, it is
$$
\HH^i(\Omega^\bullet_{X|k}) = h^i(\Omega^\bullet_{X|k}(X))
$$
where $h^i(L^\bullet)$ is the usual $i$--th cohomology of $L^\bullet$.
This holds, because $H^j(X, \Omega^p_{X|k}) = 0$ for $j > 0$ as $\Omega^p_{X|k}$ is a quasicoherent
sheaf on an affine scheme.
Theorem
Let $U = X - Y$ with $X = \AFF^n_k$ and $Y = V(f_1,\ldots,f_r)$, a closed
subscheme of $X$. It is then possible to compute
$\Hdr^i(U)$ constructively.
This is done with $D$--modules, as described in
Algorithmic Computation of de Rham Cohomology of Complements of Complex Affine Varieties,
U. Walther, J. Symb. Comput. 29, No. 4--5, 795--839 
From
On the De Rham cohomology of algebraic varieties,
R. Hartshorne, Publ. Math., Inst. Hautes Étud. Sci. 45, 5--99
we take, in the situation of the Theorem and with the additional assumption of $Y$ being smooth,
the sequence
$$
\cdots H_q(Y) \to H_q(X) \to H_q(X - Y) \to H_{q-1}(Y) \to H_{q-1}(X) \to \cdots
$$
where $H_q(W)$ stands for the de Rham homology introduced there, and the relations
(Proposition 3.4)
\begin{align*}
H_q(X) & = \Hdr^{2n - q}(X) \\
H_q(Y) & = \Hdr^{2r - q}(Y)
\end{align*}
with $\dim X = n$ and $\dim Y = r$ hold.
Together this gives
$$
\cdots \to \Hdr^{2 r - q}(Y) \to \Hdr^{2 n - q}(X) \to \Hdr^{2 n - q}(U) \to 
\Hdr^{2r - q + 1}(Y) \to \Hdr^{2 n - q + 1}(X) \to \cdots
$$
or with $s = \codim(Y, X) = n-r$:
$$
\cdots \to \Hdr^{q - 2 s}(Y) \to \Hdr^q(X) \to \Hdr^q(U) \to \Hdr^{q + 1 - 2 s}(Y) \to
\Hdr^{q+1}(X) \to \cdots
$$
But now it is $\Hdr^q(X) = \Hdr^q(\AFF^n_k) = 0$ for $q > 0$, therefore
$$
\Hdr^i(Y) = \Hdr^{i + 2 s - 1}(U)
$$
Because $\Hdr^q(U)$ by the Theorem above can be computed for all $q$, this also holds true for all
$\Hdr^i(Y)$.
It is now time to discuss the case of a general quasiprojective $X \subseteq \PSP^n_k$.
We first assume, that $X$ is smooth.
We can write
$$
X = \bar{X} \cap (\PSP^n_k - Z) = \bar{X} \cap U
$$
with closed schemes $\bar{X}, Z \subseteq \PSP^n_k$ and $Z = V(g_1,\ldots,g_s)$
where $g_i$ are homogeneous polynomials in the coordinate ring of $\PSP^n_k$.
Taking the open sets
$$
U_i = D_+(g_i) \subseteq \PSP^n_k
$$
the schemes $V_i = U_i \cap X$ are an open cover of $X$ with affine schemes, in $X$ open.
It is even $V_{i_0 \cdots i_q} = U_{i_0 \cdots i_q} \cap X \subseteq U_{i_0 \cdots i_q}$
a closed, smooth, affine subscheme of the affine scheme
$$
U_{i_0 \cdots i_q} = U_{i_0} \cap \cdots \cap U_{i_q} = D_+(g_{i_0} \cdot \cdots \cdot g_{i_q}).
$$
With a Veronese--embedding $v_d: \PSP^n_k \to \PSP^N_k$ where
$d = \prod_{\nu = 0}^q \deg g_{i_\nu}$, the map 
$$
v_d:v_d^{-1}(D_+(w)) = U_{i_0 \cdots i_q} \to D_+(w) = \AFF^N_k
$$
is a closed immersion.
In the above $w$ is a linear form in the coordinate ring of $\PSP^N_k$
derived from $\prod_{\nu=1}^q g_{i_\nu}$.
As $i: V_{i_0 \cdots i_q} \to U_{i_0 \cdots i_q}$ is a closed immersion,
the map
$$
v'_d = v_d \circ i:V_{i_0 \cdots i_q} \to U_{i_0 \cdots i_q} \xrightarrow{v_d} D_+(w) = \AFF^N_k
$$
is a closed immersion too. 
So $V_{i_0 \cdots i_q} \cong v'_d(V_{i_0 \cdots i_q})$
is even a closed, smooth subscheme of $\AFF^N_k$.
By the considerations further above $\Hdr^p(V_{i_0\cdots i_q})$ is explicitly computable.
Consider now the Čech-de Rham spectral sequence of the cover $\cover{V} = (V_i)_{i\in I}$
of $X$:
$$
E^{pq}_0 = \Cc^q(\cover{V},\Omega_{X|k}^p) = 
\prod_{i_0 < \cdots < i_q} \Gamma(V_{i_0,\ldots,i_q}, \Omega_{X|k}^p)
$$
By general theorem $E^{pq}_r$ abuts against the $\Hdr^\bullet(X)$:
$$
E^{pq}_r \Rightarrow \Hdr^\bullet(X)
$$
See for this
On the de Rham cohomology of algebraic varieties, M. Stevenson
section 4.1, p. 7.
If one forms $E^{pq}_1$ taking cohomology along the $p$--axis, the result is
$$
E^{pq}_1 = \prod_{i_0 < \cdots < i_q} \Hdr^p(V_{i_0, \ldots, i_q})
$$
Keep in mind for this that because $V_{i_0,\ldots,i_q}$ is affine, the equality
$$
\Hdr^p(V_{i_0,\ldots,i_q}) = \HH^p(\Omega^\bullet_{V_{i_0,\ldots,i_q}|k}) = 
h^p(\Gamma(V_{i_0,\ldots,i_q}, \Omega^\bullet_{V_{i_0,\ldots,i_q}|k}))
$$
holds.
As remarked above $\Hdr^p(V_{i_0,\ldots,i_q})$ is explicitly computable and, at least
theoretically, the abutment $\Hdr^\bullet(X)$ of the spectral sequences allows to be
effectively computed.
Remark
The advantage of computing with $E^{pq}_1 = \prod_{i_0 < \cdots < i_q} \Hdr^p(V_{i_0,\ldots,i_q})$ 
instead of starting with $E^{pq}_0$ directly, is, that the elements of 
$E^{pq}_0$ are infinite, but the $E^{pq}_1$ are finite $k$--modules.
The following Theorem is from Hartshorne's articles cited above (Theorem 4.4):
Theorem
Let $f:X' \to X$ be a proper map of schemes, $Y \subseteq X$ a subscheme of
$X$, and $Y' = f^{-1}(Y)$ the preimage-scheme of $Y$.
It shall hold


*

*The morphism $f$ maps $X' - Y'$ isomorphically to $X - Y$.

*Furthermore there are closed immersions $X' \to Z'$ and $X \to Z$ into smooth schemes
$Z'$, $Z$, together with a proper morphism $g:Z' \to Z$ with 
$$
\begin{CD}
X' @>>>  Z'  \\
@VfVV        @VgVV \\
X @>>>  Z
\end{CD}
$$
so that $Z' - g^{-1}(Y)$ is mapped isomorphically to $Z - Y$.
Then there is a long exact sequence in de Rham--cohomology
\begin{equation}
\cdots \to \Hdr^q(X) \to \Hdr^q(X') \oplus \Hdr^q(Y) \to \Hdr^q(Y') \to 
\Hdr^{q+1}(X) \to \cdots
\end{equation}
(end of theorem)
We can use this theorem, to compute for an arbitrary quasi-projective scheme $X \subseteq \PSP^n_k = Z$
the cohomologies $\Hdr^q(X)$.
We write $X = X_0 \cap U$ with a closed $X_0 \subseteq \PSP^n_k$ and an
open $U \subseteq \PSP^n_k$.
We find for $X_0$ a desingularization $f:X_0' \to X_0$ together with a proper
morphism $g:Z' \to Z$ with
$$
\begin{CD}
X_0' @>>>  Z'  \\
@VfVV        @VgVV \\
X_0 @>>>  Z
\end{CD}
$$
where the horizontal maps are closed immersions and $X_0'$ and $Z'$ are smooth schemes.
(See the remark after the theorem cited above in Hartshorne's article). The map
$g$ is the iterative blow-up in nonsingular subschemes.
The morphism $f:X_0' \to X_0$ fulfills the conditions of the above theorem with $X = X_0$,
$X' = X_0'$ and $Y = Y_0 \subseteq X_0$, suitable subscheme with $\dim Y_0 < \dim X_0$, 
and $Y' = Y_0' = f^{-1}(Y_0)$. One can choose $Y_0$ to be the subscheme of singular points of $X_0$, 
("strong desingularization").
Forming the base-extension of $g$ with $U \to Z$ and those of $f$ with $U \cap X_0 \to X_0$,
we have a compatible system
\begin{equation*}
\begin{CD}
X' = f^{-1}(U \cap X_0) @>>> g^{-1}(U)   \\
@Vf'VV        @VgVV \\
X @>>>  U
\end{CD}
\end{equation*}
so that the morphism $f':f^{-1}(U \cap X_0) \to X$ fulfills the conditions of the theorem above with
\begin{align}
X = X, &&  X' = f^{-1}(U \cap X_0), && Y = Y_0 \cap U, && Y' = f'^{-1}(Y_0 \cap U).
\end{align}
But $X'$ as an open part of $X_0'$ is a smooth quasi--projective scheme,
so that $\Hdr^q(X')$ is computable. Also $\Hdr^q(Y)$ and $\Hdr^q(Y')$ are computable by induction over
the dimension of the variety, for which cohomology is to be computed.
The long exact sequence from the theorem above gives therefore conditions, from which
$\Hdr^q(X)$ can be computed:
\begin{multline}
\Hdr^{q}(X') \oplus \Hdr^q(Y) \to \Hdr^q(Y') \to \Hdr^{q+1}(X) \to \\
\to \Hdr^{q+1}(X') \oplus \Hdr^{q+1}(Y) \to \Hdr^{q+1}(Y')
\end{multline}
So we can compute $\Hdr^q(X)$ for all quasiprojective 
$X \subseteq \PSP^n_k$.
A: Note that the question was not about the cohomology, but about the spaces of differential forms.
Similarly to the given recipe for $\Omega^1(X)$, the space of differential forms of higher order on an affine variety $X=\operatorname{Spec}[x_1,\ldots,x_n]/(f_1,\ldots,f_M)$ is the quotient module
$$
\Omega^k(X) = \frac{\bigoplus_{|I|=k} O(X) d x_I}{\left(\sum_{i=1}^n \frac{\partial f_m}{\partial x_i} d x_i \wedge d x_J \Big| |J|=k-1,\, m=1,\ldots,M\right)}.
$$
Here $I$, $J$ are subsets of $\{1,\ldots,n\}$ and for such a subset $I=\{i_1,\ldots,i_k\}$ with $i_1<\cdots<i_k$ we denote $d x_I = d x_{i_1} \wedge \cdots \wedge d x_{i_k}$.
This is simply by unwrapping the definition $\Omega^k = \wedge^k \Omega^1$.
