CoVaR/dCoVaR modelling using bivariate DCC-GJR-GARCH

Question

For the several weeks, I have been looking for a way to calculate and display the results of my DCC-GJR-GARCH model to picture a dynamic relationship between daily return of, let's say for example, two assets. I read many reference papers and they all got similar procedure, it's something like this, (1) first step is to estimating a univariate GARCH model of each assets, (2) second is to estimate a bivariate GARCH with DCC specification, and (3) third is to numerically solve double integral method to obtain CoVaR/dCoVaR (https://www.sciencedirect.com/science/article/abs/pii/S0378426613001155). Now I lost it at the third step. I'm using R and this is the code that I have been using.

### Set and load libraries

library(tidyverse)
library(quantmod)
library(rugarch)
library(rmgarch)
library(ggplot2)

### Download dataset
### In this case I'm using Bitcoin and SP500 to check spillover effect from Bitcoin

d1=getSymbols("BTC-USD",auto.assign=F)|>Ad()|>ROC()|>fortify.zoo()
d2=getSymbols("^GSPC",auto.assign=F)|>Ad()|>ROC()|>fortify.zoo()
colnames(d1)=c("Index","BTC")
colnames(d2)=c("Index","SP500")
d3=na.omit(merge(d1,d2,all=F))

### Step 1 (univariate GJR-GARCH model)

m1=ugarchspec(variance.model=list(model="gjrGARCH",garchOrder=c(1,1)),
              mean.model=list(armaOrder=c(0,0),include.mean=T),
              distribution.model="std")
m2=ugarchfit(spec=m1,data=d3$BTC)
m3=ugarchfit(spec=m1,data=d3$SP500)

d4=fitdist(distribution='std',x=d3$BTC)$pars[3]
d5=qdist(distribution='std',shape=d4,p=0.05)
d6=fitdist(distribution='std',x=d3$SP500)$pars[3]
d7=qdist(distribution='std',shape=d6,p=0.05)

d8=m2@fit$sigma*d5
d9=m3@fit$sigma*d7
d10=cbind(d3,d8,d9)
colnames(d10)[4]="BTCvol"
colnames(d10)[5]="SP500vol"

f1=ggplot(data=d10,aes(x=Index))+geom_line(aes(y=BTCvol,color="BTC"))+
  geom_line(aes(y=SP500vol,color="SP500"))+ylab("")+xlab("")+
  scale_color_manual(name="Volatility Series", values=c("BTC"="red","SP500"="darkblue"))+
  theme_light()+theme(legend.position="top")

### Step 2 (DCC-GARCH specification)

d11=d10[,c("BTC","SP500")]
m4=dccspec(uspec=multispec(replicate(2,m1)),dccOrder=c(1,1),distribution="mvt")
m5=dccfit(m4,data=d11)
d12=m5@model$sigma

That was the last line of code that I wrote, I'm stuck as I don't know what else to do. While my goal is to make something like this:

And also this:

That was from a paper by Zulu (https://ujcontent.uj.ac.za/esploro/outputs/graduate/The-contribution-of-South-African-Insurance/9917608607691). The point is, my main goal is to examine how much risk cryptocurrencies contribute to a stock market in my country, and compared it with other nation's stock market, And make some charts like the example above. Also I don't quite understand the term "system" mentioned in papers above, conceptually I understand, but technically I don't quite get how to measure the "system"'s CoVaR.